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

[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 ?