Electro-hydraulic actuator system simulation generates complex result under open-loop simulation mode

I have generated simulation datas of electro-hydraulic actuator (EHA) using the explicit physical model based on MATLAB's ode solver. After that, I want to simulate the EHA based on discrete time model which is obtained from the first-order euler discretization of the continue time countpart . However, the simulation result of the discrete time model generates complex result caused by incorrect negative pressure. How to resolve this issue, I suspect that the error comes from the discretization method.

% Construct a struct for storing EHA parameters
EHAParams = struct;
EHAParams.A1 = 1.96e-3; % Hydraulic cylinder piston side area
EHAParams.A2 = 1.65e-3; % Hydraulic cylinder rod side area
EHAParams.L = 0.5; % Cylinder stroke
EHAParams.V10 = 5e-5; % Initial Volume V10
EHAParams.V20 = 5e-5; % Initial Volume V20

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

EHAParams.rho = 850; % Fluid density
EHAParams.Cd = 0.6; % Discharge coefficient
EHAParams.Wv = 20e-3; % Servo valve area gradient
EHAParams.ku = 5e-5; % Spool displacement gain
EHAParams.Cil = 1e-8; % hydraulic cylinder internal leakage coeff
EHAParams.Cel = 1e-8; % hydraulic cylinder external leakage coeff

% PID controller parameters
EHAParams.PID.Kp = 300;
EHAParams.PID.Ki = 200;
EHAParams.PID.Kd = 20; % 

% reference trajectory and velocity
EHAParams.refPos = @(t) 0.1 + 0.05*sin(pi*t);
EHAParams.refVel = @(t) 0.05*pi*cos(pi*t);

% Consistent initial conditions
P10 = 1e7;
P20 = 1e7;

X0 = [0.1;0;P10;P20;0]; % initial state of EHA

startTime = 0;
endTime = 5;

F = ode("ODEFcn",@(t,x,p) EHAODEs(t,x,p),"InitialTime",startTime,...
        "InitialValue",X0,"Parameters",EHAParams,"Solver","ode15s");
F.RelativeTolerance = 1e-6;
solFcn = solutionFcn(F,startTime,endTime);

% Synthetic data generation using odes
sampleFreq = 1000; % sample frequency 
sampleTime = startTime:1/sampleFreq:endTime;
trueStates = solFcn(sampleTime);

extInput = EHAParams.PID.Kp*(EHAParams.refPos(sampleTime) - trueStates(1,:)) + ...
           EHAParams.PID.Kd*(EHAParams.refVel(sampleTime) - trueStates(2,:)) + ...
           EHAParams.PID.Ki*trueStates(end,:);
plot(sampleTime,trueStates(1,:),sampleTime,EHAParams.refPos(sampleTime),'--');
legend(["Output position","Reference position"]);

plot(sampleTime,extInput,"LineWidth",2);

predStates = zeros(size(trueStates(1:4,:)));
predStates(:,1) = trueStates(1:4,1);

for p = 2:size(predStates,2)
    predStates(:,p) = stateTransitionEHA(predStates(:,p-1),1/sampleFreq,extInput(p-1),EHAParams);    
end

isreal(predStates)

ans = logical
   0

% discrete time counterpart of EHAODEs
function Xnext = stateTransitionEHA(Xcur,dt,Ucur,P)
    % discrete time state transition model of EHA
    %
    % Input arguments: Xcur -> current state of EHA
    %                  dt -> discrete time step
    %                  Ucur -> current control input of EHA
    %                  P -> parameters struct of EHA system

    % Output arguments: Xnext -> predicted state of EHA
    % Note: state of EHA contains [position;velocity;pressure-1;pressure-2]

    Xv = P.ku * Ucur; % displacement of spool

    dX = zeros(4,1);
    dX(1) = Xcur(2);
    dX(2) = (P.A1*Xcur(3) - P.A2*Xcur(4) - ...
            P.k*Xcur(1) - P.b*Xcur(2) - ...
            P.FL) / P.m;

    if Xv >= 0
        Q1 = P.Cd * P.Wv * Xv * sqrt(2/P.rho * (P.Ps - Xcur(3)));
        Q2 = P.Cd * P.Wv * Xv * sqrt(2/P.rho * Xcur(4));

    else
        Q1 = P.Cd * P.Wv * Xv * sqrt(2/P.rho * Xcur(3));
        Q2 = P.Cd * P.Wv * Xv * sqrt(2/P.rho * P.Ps - Xcur(4));

    end

    V1 = P.V10 + P.A1 * Xcur(1); 
    V2 = P.V20 + P.A2 * (P.L - Xcur(1));

    dX(3) = (P.beta / V1) * (Q1 - P.A1 * Xcur(2) - P.Cil * (Xcur(3)-Xcur(4)));
    dX(4) = (P.beta / V2) * (-Q2 + P.A2 * Xcur(2) + P.Cil * (Xcur(3) - Xcur(4)) - P.Cel * Xcur(4));

    Xnext = Xcur + dt * dX;
end

% state space model of EHA
function dX = EHAODEs(t,X,P)
    % ODEs of EHA
    % Input arguments: t -> current time
    %                  X -> current state
    %                  P -> parameter struct of EHA

    % Output arguments: dX -> differential of state X
    % X == [pos;vel;P1;P2;errInt] under PID closed-loop control

    dX = zeros(5,1);

    u = P.PID.Kp * (P.refPos(t) - X(1)) + ...
        P.PID.Kd * (P.refVel(t) - X(2)) + ...
        P.PID.Ki * X(5);
    xv = P.ku * u; % servo valve spool diaplacement

    dX(1) = X(2);
    dX(5) = P.refPos(t) - X(1); % position error

    dX(2) = (P.A1*X(3) - P.A2*X(4) - ...
            P.k*X(1) - P.b*X(2) - ...
            P.FL - sin(pi*t)) / P.m; % mechanical dynamic equation

    % differential equations of EHA pressures

    if xv >= 0
        Q1 = P.Cd * P.Wv * xv * sqrt(2/P.rho * (P.Ps - X(3)));
        Q2 = P.Cd * P.Wv * xv * sqrt(2/P.rho * X(4));
    else
        Q1 = P.Cd * P.Wv * xv * sqrt(2/P.rho * X(3));
        Q2 = P.Cd * P.Wv * xv * sqrt(2/P.rho * (P.Ps - X(4)));
    end
    
    V1 = P.V10 + P.A1 * X(1); 
    V2 = P.V20 + P.A2 * (P.L - X(1));

    dX(3) = (P.beta / V1) * (Q1 - P.A1 * X(2) - P.Cil * (X(3)-X(4)));
    dX(4) = (P.beta / V2) * (-Q2 + P.A2 * X(2) + P.Cil * (X(3) - X(4)) - P.Cel * X(4));
end

Answer 1

[Chuguang Pan]

I have find that the complex result issue appeared in the single step forward prediction is induced by the discretization error of continuous state space model. Since the first order Euler forward discretization is not appropriate for the nonlinear EHA system, it is prone to generate cumulative error in the continuous single step forward prediction. After using higher precision discretization method, the forward prediction result is consistent to the real state.

Answer 2

[Pat Gipper]

I believe the hydraulic cylinder external leakage coefficient is too large and is pulling the cylinder pressures down near the tank pressure (0). Try reducing it by a factor of 10^4.

[Chuguang Pan]

@Pat Gipper. Thanks for your answer. After setting the external leakage coefficient from 1e-8 to 1e-12, the discrete time one-step forward prediction still generating complex result caused by negative predicted pressure.

[Pat Gipper]

Hi @Chuguang Pan. Sorry my first suggestion didn't help. I ran into a couple issues when I tried to use your original script files. Eventually I set the external leakage coefficient to zero. Once I did this the solver produced an error very early in the solution when the time step became too small. Making the cylinder areas balanced by setting EHAParams.A2 = EHAParams.A1 fixed this problem and allowed the solution to run to completion.

[Pat Gipper]

I think that when adding valve leakage terms to the cylinders volumes that you should try to keep them balanced, meaning that if you have leakage to the tank then you should have an equivalent amount of leakage from the supply. This should be on both the extend and retract sides of the cylinder. The whole point of this is to keep the cylinder neutral pressures balanced near 1/2 the supply pressure.

Before exploring the closed loop system response you should determine that cylinder neutral pressures will settle at desired values. Do this by setting the reference commands to a fixed position and zero velocity and also remove all external loads. Once you are happy with the neutral pressures then go ahead and explore system performance.

[Chuguang Pan]

@Pat Gipper. I am confused about how to determine the cylinder neutral pressure by settling the reference commands to a fixed position and zero velocity, can you give me some concrete code examples through revising my code ?

[Pat Gipper]

@Chuguang Pan I have attached a modification to your code that demonstrates what I mean. Lines 36 through 44 modifies the parameters to remove external loads and command a fixed reference position. I also needed the modfication shown on line 104 to remove a sinusoidal load being applied in the differential equations. I was able to make this work with the unbalanced area shown in line 38. Smaller values cause a problem with the solution time step. The plot below shows the steady-state pressures in the cylinder. As you can see the have settled at about 1/2 the supply pressure, differing because of the unbalanced area. The steady states values are what are called the neutral pressure.

[Chuguang Pan]

@Pat Gipper. Thanks for your explanation. I have find that the open-loop one-step forward prediction generates significant deviation from the trueStates when I use discrete time state transition model stateTransitionEHA. Do you konw how to resolve this accumulated one-step prediction error.