Electro-hydraulic actuator system simulating based on ode object generates error

I want to simulate the electro-hydraulic actuator (EHA) system using MATLAB's ode object. The schematic diagram of the EHA is shown as follows.

However, the code does not generate effective result. I think the error is induced by the stiffness of the EHA system, what should I do to obtain the effective result?

% EHA system parameters initialization
EHAParams = struct;
EHAParams.A1 = 1.96e-3; % Hydraulic cylinder piston side area
EHAParams.A2 = 1.65e-4; % Hydraulic cylinder rod side area
EHAParams.L = 0.5; % Cylinder stroke
EHAParams.V10 = EHAParams.A1 * EHAParams.L / 2; % Initial Volume V10 [mid position]
EHAParams.V20 = EHAParams.A2 * EHAParams.L / 2; % Initial Volume V20 [mid position]

EHAParams.beta = 1700e6; % Effective bulk modulus
EHAParams.b = 100; % Friction coefficient
EHAParams.k = 20; % Stiffness coefficient
EHAParams.m = 50; % Load mass
EHAParams.Pt = 0; % Tank pressure
EHAParams.Ps = 21e6; % Supply pressure
EHAParams.FL = 100; % Constant load

EHAParams.rho = 850; % Fluid density
EHAParams.Cd = 0.62; % Discharge coefficient
EHAParams.Wv = pi*1e-3/4; % Servo valve area gradient
EHAParams.ku = 2e-4; % Spool displacement gain

% PID controller parameters
 EHAParams.PID.Kp = 200;
 EHAParams.PID.Ki = 50;
 EHAParams.PID.Kd = 1; 
 % reference position and velocity cmd for PID
 EHAParams.refPos = @(t) 0.1 + 0.05*sin(pi*t);
 EHAParams.refVel = @(t) 0.05*pi*cos(pi*t);

% Simulate the EHA system using ode object
X0 = [0;0;1e7;1e7;0]; % initial state [pos;vel;P1;P2;error_int]
startTime = 0;
endTime = 10;

F = ode(ODEFcn=@(t,x,p) stateTransitionEHA(t,x,p),InitialTime=startTime,InitialValue=X0,...
        Parameters=EHAParams,Solver="ode23t");
F.SolverOptions.OutputFcn = @odeplot;
F.SolverOptions.OutputSelection = [3 4];

sol = solve(F,startTime,endTime);

Error using odePlotImpl (line 61)
Error updating the ODEPLOT window. Solution data may have been corrupted. Argument Y cannot be complex.

Error in odeplot (line 40)
status = odePlotImpl(fun,flag,"odeplot",SetAxis=setAxis);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in ode23t (line 758)
stop = feval(outputFcn,tout_new,yout_new(outputs,:),'',outputArgs{:});
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in matlab.ode.internal.LegacySolver/runSolver (line 672)
[varargout{1:nout}] = obj.SolverFcn(obj.ODEFcn, tspan, ...
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in matlab.ode.internal.LegacySolver/solve (line 141)
[t,y] = obj.runSolver(obj.InitialTime,t2);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in ode/solveBetween (line 1009)
[tR,yR] = solver.solve(t0,tmax,pps);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in ode/solve (line 388)
sol = obj.solveBetween(t1,t2,refine);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

function dX = stateTransitionEHA(t,X,Params)
    % State transition model of EHA 
    % Input arguments: t->current time
    %                  X->current state                 
    %                  Params -> struct contains EHA parameters
    % Output arguments: dX -> state differential

     H = @(u) max(0,sign(u)); % Heaviside function H
     
     dX = zeros(5,1); % Preallocation state derivative

     dX(1) = X(2);
     dX(2) = (-Params.k*X(1) - Params.b*X(2) + Params.A1*X(3) - Params.A2*X(4) - ...
               Params.FL)/Params.m; % 
     u = Params.PID.Kp*(Params.refPos(t) - X(1)) + Params.PID.Kd*(Params.refVel(t) - X(2)) + ...
         Params.PID.Ki*X(5); % current control input
     % ------------------------------------------------------- %
     kq = Params.Cd * Params.Wv * sqrt(2/Params.rho); 
     g1 = H(u)*sqrt(Params.Ps - X(3)) + H(-u)*sqrt(X(3));
     g2 = H(u)*sqrt(X(4)) + H(-u)*sqrt(Params.Ps - X(4));
     h1 = 1 / (Params.V10 + Params.A1 * X(1));
     h2 = 1 / (Params.V20 - Params.A2 * X(1));

     dX(3) = Params.beta * h1 * (kq * Params.ku * g1 * u - Params.A1 * X(2) - ...
                                 2e-6 * (X(3) - X(4)));
     dX(4) = Params.beta * h2 * (-kq * Params.ku * g2 * u + Params.A2 * X(2) + ...
                                 2e-6 * (X(3) - X(4)));
     dX(5) = Params.refPos(t) - X(1); % tracking error (ref - pos)
end

Answer 1

[Satyam]

Hi,

The error occurs because during integration some terms inside sqrt can become negative, which produces complex values and causes odeplot to fail. Constraining the square-root arguments and preventing singular cylinder volumes allows the solver to run properly. A stiff solver is also more appropriate for this hydraulic system. An example corrected implementation is shown below.

% EHA system parameters initialization

EHAParams = struct;

EHAParams.A1 = 1.96e-3;

EHAParams.A2 = 1.65e-4;

EHAParams.L = 0.5;

EHAParams.V10 = EHAParams.A1 * EHAParams.L / 2;

EHAParams.V20 = EHAParams.A2 * EHAParams.L / 2;

EHAParams.beta = 1700e6;

EHAParams.b = 100;

EHAParams.k = 20;

EHAParams.m = 50;

EHAParams.Pt = 0;

EHAParams.Ps = 21e6;

EHAParams.FL = 100;

EHAParams.rho = 850;

EHAParams.Cd = 0.62;

EHAParams.Wv = pi*1e-3/4;

EHAParams.ku = 2e-4;

% PID controller parameters

EHAParams.PID.Kp = 200;

EHAParams.PID.Ki = 50;

EHAParams.PID.Kd = 1;

% reference signals

EHAParams.refPos = @(t) 0.1 + 0.05*sin(pi*t);

EHAParams.refVel = @(t) 0.05*pi*cos(pi*t);

% Initial state [pos; vel; P1; P2; error_int]

X0 = [0;0;1e7;1e7;0];

startTime = 0;

endTime = 10;

% Use stiff solver

F = ode(ODEFcn=@(t,x,p) stateTransitionEHA(t,x,p), ...

InitialTime=startTime, ...

InitialValue=X0, ...

Parameters=EHAParams, ...

Solver="ode15s");

F.SolverOptions.OutputFcn = @odeplot;

F.SolverOptions.OutputSelection = [3 4];

sol = solve(F,startTime,endTime);

function dX = stateTransitionEHA(t,X,Params)

H = @(u) max(0,sign(u));

dX = zeros(5,1);

% Mechanical dynamics

dX(1) = X(2);

dX(2) = (-Params.k*X(1) ...

-Params.b*X(2) ...

+Params.A1*X(3) ...

-Params.A2*X(4) ...

-Params.FL)/Params.m;

% PID controller

u = Params.PID.Kp*(Params.refPos(t)-X(1)) ...

+ Params.PID.Kd*(Params.refVel(t)-X(2)) ...

+ Params.PID.Ki*X(5);

% Flow equations

kq = Params.Cd*Params.Wv*sqrt(2/Params.rho);

g1 = H(u)*sqrt(max(Params.Ps-X(3),0)) ...

+ H(-u)*sqrt(max(X(3),0));

g2 = H(u)*sqrt(max(X(4),0)) ...

+ H(-u)*sqrt(max(Params.Ps-X(4),0));

% Prevent volume singularities

V1 = max(Params.V10 + Params.A1*X(1),1e-9);

V2 = max(Params.V20 - Params.A2*X(1),1e-9);

h1 = 1/V1;

h2 = 1/V2;

% Pressure dynamics

dX(3) = Params.beta*h1*(kq*Params.ku*g1*u ...

- Params.A1*X(2) ...

- 2e-6*(X(3)-X(4)));

dX(4) = Params.beta*h2*(-kq*Params.ku*g2*u ...

+ Params.A2*X(2) ...

+ 2e-6*(X(3)-X(4)));

% Integral error

dX(5) = Params.refPos(t) - X(1);

end

Since electro-hydraulic systems are typically stiff, use a stiff solver such as ode15s or ode23s instead of ode23t, and optionally tighten tolerances through F.SolverOptions. MATLAB documentation on stiff solvers explains when these methods are more appropriate: https://www.mathworks.com/help/matlab/math/solve-stiff-odes.html.

I hope it resolves your query.

[Chuguang Pan]

@Satyam. Thanks for your answer. However, the simulation results are not resonable since the hydraulic cylinder's output displacement is very large.

% EHA system parameters initialization
EHAParams = struct;
EHAParams.A1 = 1.96e-3;
EHAParams.A2 = 1.65e-4;
EHAParams.L  = 0.5;
EHAParams.V10 = EHAParams.A1 * EHAParams.L / 2;
EHAParams.V20 = EHAParams.A2 * EHAParams.L / 2;
EHAParams.beta = 1700e6;
EHAParams.b = 100;
EHAParams.k = 20;
EHAParams.m = 50;
EHAParams.Pt = 0;
EHAParams.Ps = 21e6;
EHAParams.FL = 100;
EHAParams.rho = 850;
EHAParams.Cd = 0.62;
EHAParams.Wv = pi*1e-3/4;
EHAParams.ku = 2e-4;
% PID controller parameters
EHAParams.PID.Kp = 200;
EHAParams.PID.Ki = 50;
EHAParams.PID.Kd = 1;
% reference signals
EHAParams.refPos = @(t) 0.1 + 0.05*sin(pi*t);
EHAParams.refVel = @(t) 0.05*pi*cos(pi*t);

% Initial state [pos; vel; P1; P2; error_int]
X0 = [0;0;1e7;1e7;0];
startTime = 0;
endTime = 10;

% Use stiff solver
F = ode(ODEFcn=@(t,x,p) stateTransitionEHA(t,x,p), ...
        InitialTime=startTime, ...
        InitialValue=X0, ...
        Parameters=EHAParams, ...
        Solver="ode15s");

F.SolverOptions.OutputFcn = @odeplot;
F.SolverOptions.OutputSelection = 1;

sol = solve(F,startTime,endTime);

function dX = stateTransitionEHA(t,X,Params)
H = @(u) max(0,sign(u));
dX = zeros(5,1);
% Mechanical dynamics
dX(1) = X(2);
dX(2) = (-Params.k*X(1) ...
        -Params.b*X(2) ...
        +Params.A1*X(3) ...
        -Params.A2*X(4) ...
        -Params.FL)/Params.m;
% PID controller
u = Params.PID.Kp*(Params.refPos(t)-X(1)) ...
  + Params.PID.Kd*(Params.refVel(t)-X(2)) ...
  + Params.PID.Ki*X(5);
% Flow equations
kq = Params.Cd*Params.Wv*sqrt(2/Params.rho);
g1 = H(u)*sqrt(max(Params.Ps-X(3),0)) ...
   + H(-u)*sqrt(max(X(3),0));
g2 = H(u)*sqrt(max(X(4),0)) ...
   + H(-u)*sqrt(max(Params.Ps-X(4),0));
% Prevent volume singularities
V1 = max(Params.V10 + Params.A1*X(1),1e-9);
V2 = max(Params.V20 - Params.A2*X(1),1e-9);
h1 = 1/V1;
h2 = 1/V2;
% Pressure dynamics
dX(3) = Params.beta*h1*(kq*Params.ku*g1*u ...
        - Params.A1*X(2) ...
        - 2e-6*(X(3)-X(4)));
dX(4) = Params.beta*h2*(-kq*Params.ku*g2*u ...
        + Params.A2*X(2) ...
        + 2e-6*(X(3)-X(4)));
% Integral error
dX(5) = Params.refPos(t) - X(1);
end