[Math] Why is the state space model of a hydraulic servo system unstable

Question

I have those equations:

Flow valve equation – in:
$$Q_1 = C_qax\sqrt{\frac{2}{\rho}(P_s -P_1)}$$
Flow valve equation – Out:
$$Q_2 = C_qax\sqrt{\frac{2}{\rho}(P_2)}$$
Hydraulic actuator equation – flow in:
$$Q_1 = A_1\dot{y} + \frac{V_1}{\beta}\dot{P_1}$$
Hydraulic actuator equation – flow Out:
$$Q_2 = A_2\dot{y} + \frac{V_2}{\beta}\dot{P_2}$$
Hydraulic actuator equation – Force:
$$M\ddot{y} = A_1P_1 – A_2P_2 – C_f\dot{y} – F_L$$

I want to combine them all into a state space model. If you wonder what they are, they are hydraulic equations for a servo actuator.

Where:

• $M$ is mass $\text{kg}$
• $y$ is position of the mass
• $A$ is area $\text{m^2}$
• $P$ is pressure $\text{Pa}$
• $F_L$ is load $\text{M}$
• $C_f$ is friction $\text{Ns/m}$
• $C_q$ is flow coffecient for the valve, most of the time $C_q = 0.67$
• $a$ is area gradient for the valve
• $x$ is position of valve $\text{m}$
• $\beta$ is bulk modulus of the oil $\text{Pa}$
• $V$ is the volume of the chamber of the actuator $\text{m^3}$

Anyway! I want to combine them into a state space model. So I need to linearize the valve flow equations:

$$Q_1 = C_qax\sqrt{\frac{2}{\rho}(P_s -P_1)}$$
$$Q_2 = C_qax\sqrt{\frac{2}{\rho}(P_2)}$$

Those will be:
Note that this will be negative!

$$p_1 = \frac{\partial Q_1}{\partial P_1} = -\frac{C_qax}{2\sqrt{\rho(P_s – P_1)}}$$

$$q_1 = \frac{\partial Q_1}{\partial x} = C_qa\sqrt{\frac{2}{\rho}(P_s -P_1)}$$

and

$$p_2 = \frac{\partial Q_2}{\partial P_2} = \frac{C_qax}{2\sqrt{\rho(P_2)}}$$

$$q_2 = \frac{\partial Q_2}{\partial x} = C_qa\sqrt{\frac{2}{\rho}(P_2)}$$

So our linearized valve flow equations are then:

$$Q_1 = q_1x – p_1P_1$$
$$Q_2 = q_2x + p_2P_2$$

And now we will combine:

$$ q_1x – p_1P_1 = A_1\dot{y} + \frac{V_1}{\beta}\dot{P_1}$$
$$ q_2x + p_2P_2 = A_2\dot{y} + \frac{V_2}{\beta}\dot{P_2}$$
$$M\ddot{y} = A_1P_1 – A_2P_2 – C_f\dot{y} – F_L$$

The purpose to do a state space is that we need to write the system on the first order equations:

So one method is to write
$$y = y_1, \dot{y_1} = y_2 = \dot{y}, \ddot{y} = \dot{y_2}$$

The equations will look like:

$$ q_1x – p_1P_1 = A_1y_2 + \frac{V_1}{\beta}\dot{P_1}$$
$$ q_2x + p_2P_2 = A_2y_2 + \frac{V_2}{\beta}\dot{P_2}$$
$$\dot{y_1} = y_2$$
$$M\dot{y_2} = A_1P_1 – A_2P_2 – C_fy_2 – F_L$$

And we move around the equations:
$$ \dot{P_1} = \frac{\beta q_1x}{V_1} – \frac{\beta p_1P_1}{V_1} – \frac{\beta A_1y_2}{V_1}$$

$$ \dot{P_2} = \frac{\beta q_2x}{V_2} + \frac{\beta p_2P_2}{V_2} – \frac{\beta A_2y_2}{V_2}$$

$$\dot{y_1} = y_2$$
$$\dot{y_2} = \frac{A_1P_1}{M} – \frac{A_2P_2}{M} – \frac{C_fy_2}{M} – \frac{F_L}{M}$$

Now we create our state space model:

$$\dot{x} = Ax + Bu \\ y = Cx + Du$$

And that will be:

$$\begin{bmatrix}
\dot{P_1}\\
\dot{P_2}\\
\dot{y_1}\\
\dot{y_2}
\end{bmatrix} = \begin{bmatrix}
-\frac{\beta p_1}{V_1} & 0 & 0 & -\frac{\beta A_1}{V_1}\\
0& \frac{\beta p_2}{V_2} &0 & -\frac{\beta A_2}{V_2}\\
0 & 0 &0 &1 \\
\frac{A_1}{M}& -\frac{A_2}{M} &0 & -\frac{C_f}{M}
\end{bmatrix}\begin{bmatrix}
P_1\\
P_2\\
y_1\\
y_2
\end{bmatrix} + \begin{bmatrix}
\frac{\beta q_1}{V_1} & 0\\
\frac{\beta q_2}{V_2} & 0\\
0 & 0\\
0 & -\frac{1}{M}
\end{bmatrix}\begin{bmatrix}
x\\
F_L
\end{bmatrix} \\
\begin{bmatrix}
y
\end{bmatrix} = \begin{bmatrix}
0 &0 &1 &0
\end{bmatrix}\begin{bmatrix}
P_1\\
P_2\\
y_1\\
y_2
\end{bmatrix}+\begin{bmatrix}
0 &0
\end{bmatrix}\begin{bmatrix}
x\\
F_L
\end{bmatrix}$$

Question:

It still gives a unstable system – Why?

>> A = [-3 0 0 -5; 0 2 0 -2; 0 0 0 1; 3 -5 0 -2]
A =

  -3   0   0  -5
   0   2   0  -2
   0   0   0   1
   3  -5   0  -2

>> B = [3 0; 6 0; 0 0; 0 -5]
B =

   3   0
   6   0
   0   0
   0  -5

>> eig(A)
ans =

   0.00000 + 0.00000i
   3.30174 + 0.00000i <-- Positive = Unstable
  -3.15087 + 3.44655i
  -3.15087 - 3.44655i

>>

Edit:

I try hydraulic values.

A1 = 2.5*10^(-3);
A2 = A1;
Be = 1.0*10^9;
M = 15;
V1 = 1.0*10^(-3);
V2 = V1;
Cq = 0.67;
a = 0.05;
Cf = 10;
rho = 900;
x0 = 0.003;
Ps0 = 50*10^6;
P10 = 1*10^6;
P20 = 1*10^6;

p1 = Cq*a*x0/(2*sqrt(rho*(Ps0 -P10)))
q1 = Cq*a*sqrt(2/rho*(Ps0 - P10))

p2 = Cq*a*x0/(2*sqrt(rho*(P20)))
q2 = Cq*a*sqrt(2/rho*(P20))

A = [-Be*p1/V1 0 0 -Be*A1/V1; 
     0 Be*p2/V2 0 -Be*A2/V2; 
     0 0 0 1;
     A1/M -A2/M 0 -Cf/M]

B = [Be*q1/V1 0; 
     Be*q2/V1 0;
     0 0;
     0 -1/M]

eig(A) % Unstable pole  

>> eig(A)
ans =

   0.0000e+00 + 0.0000e+00i
   1.8759e+03 + 0.0000e+00i <-- Positive = Unstable
  -2.2045e+02 + 6.1378e+02i
  -2.2045e+02 - 6.1378e+02i

Edit2:

I have combined leakage $C_l$, which is the efficiency factor of the cylinder:

$$ q_1x – p_1P_1 = A_1\dot{y} + \frac{V_1}{\beta}\dot{P_1} + C_l(P_1 – P_2)$$
$$ q_2x + p_2P_2 = A_2\dot{y} + \frac{V_2}{\beta}\dot{P_2} + C_l(P_2 – P_1)$$
$$M\ddot{y} = A_1P_1 – A_2P_2 – C_f\dot{y} – F_L$$

And –>

$$ \dot{P_1} = \frac{\beta q_1x}{V_1} – \frac{\beta p_1P_1}{V_1} – \frac{\beta A_1y_2}{V_1} -\frac{\beta C_lP_1}{V_1} + \frac{\beta C_lP_2}{V_1}$$

$$ \dot{P_2} = \frac{\beta q_2x}{V_2} + \frac{\beta p_2P_2}{V_2} – \frac{\beta A_2y_2}{V_2} -\frac{\beta C_lP_2}{V_2} + \frac{\beta C_lP_1}{V_2}$$

And our state space model will be this. I also changed this $\frac{\beta q_2}{V_2}$ into $-\frac{\beta q_2}{V_2}$ because when we open the valve, the other side of the cylinder, should lose pressure.

$$\begin{bmatrix}
\dot{P_1}\\
\dot{P_2}\\
\dot{y_1}\\
\dot{y_2}
\end{bmatrix} = \begin{bmatrix}
-\frac{\beta (p_1 -C_l)}{V_1} & \frac{\beta C_l}{V_1} & 0 & -\frac{\beta A_1}{V_1}\\
\frac{\beta C_l}{V_2}& \frac{\beta (p_2-C_l)}{V_2} &0 & -\frac{\beta A_2}{V_2}\\
0 & 0 &0 &1 \\
\frac{A_1}{M}& -\frac{A_2}{M} &0 & -\frac{C_f}{M}
\end{bmatrix}\begin{bmatrix}
P_1\\
P_2\\
y_1\\
y_2
\end{bmatrix} + \begin{bmatrix}
\frac{\beta q_1}{V_1} & 0\\
-\frac{\beta q_2}{V_2} & 0\\
0 & 0\\
0 & -\frac{1}{M}
\end{bmatrix}\begin{bmatrix}
x\\
F_L
\end{bmatrix} \\
\begin{bmatrix}
y
\end{bmatrix} = \begin{bmatrix}
0 &0 &1 &0
\end{bmatrix}\begin{bmatrix}
P_1\\
P_2\\
y_1\\
y_2
\end{bmatrix}+\begin{bmatrix}
0 &0
\end{bmatrix}\begin{bmatrix}
x\\
F_L
\end{bmatrix}$$

Answer 1

Here is the answer:

Flow valve and cylinder equation for the cap-side are the same:

$$ q_1x - p_1P_1 = A_1\dot{y} + \frac{V_1}{\beta}\dot{P_1} + C_l(P_1 - P_2)$$

If the cylinder's pressure in the rod-side, increase, then the flow would stop. So instead of this:

$$ q_2x + p_2P_2 = A_2\dot{y} + \frac{V_2}{\beta}\dot{P_2} + C_l(P_2 - P_1)$$

It must be this: ( I changed leakage too) $$ q_2x + p_2P_2 = A_2\dot{y} - \frac{V_2}{\beta}\dot{P_2} + C_l(P_1 - P_2)$$

The force equation is the same: $$M\ddot{y} = A_1P_1 - A_2P_2 - C_f\dot{y} - F_L$$

The result:

$$ \dot{P_1} = \frac{\beta q_1x}{V_1} - \frac{\beta p_1P_1}{V_1} - \frac{\beta A_1y_2}{V_1} -\frac{\beta C_lP_1}{V_1} + \frac{\beta C_lP_2}{V_1}$$

$$ \dot{P_2} = -\frac{\beta q_2x}{V_2} - \frac{\beta p_2P_2}{V_2} + \frac{\beta A_2y_2}{V_2} -\frac{\beta C_lP_2}{V_2} + \frac{\beta C_lP_1}{V_2}$$

$$\begin{bmatrix} \dot{P_1}\\ \dot{P_2}\\ \dot{y_1}\\ \dot{y_2} \end{bmatrix} = \begin{bmatrix} -\frac{\beta (p_1 + C_l)}{V_1} & \frac{\beta C_l}{V_1} & 0 & -\frac{\beta A_1}{V_1}\\ \frac{\beta C_l}{V_2}& -\frac{\beta (p_2+C_l)}{V_2} &0 & \frac{\beta A_2}{V_2}\\ 0 & 0 &0 &1 \\ \frac{A_1}{M}& -\frac{A_2}{M} &0 & -\frac{C_f}{M} \end{bmatrix}\begin{bmatrix} P_1\\ P_2\\ y_1\\ y_2 \end{bmatrix} + \begin{bmatrix} \frac{\beta q_1}{V_1} & 0\\ -\frac{\beta q_2}{V_2} & 0\\ 0 & 0\\ 0 & -\frac{1}{M} \end{bmatrix}\begin{bmatrix} x\\ F_L \end{bmatrix} \\ \begin{bmatrix} y \end{bmatrix} = \begin{bmatrix} 0 &0 &1 &0 \end{bmatrix}\begin{bmatrix} P_1\\ P_2\\ y_1\\ y_2 \end{bmatrix}+\begin{bmatrix} 0 &0 \end{bmatrix}\begin{bmatrix} x\\ F_L \end{bmatrix}$$

So, I did a simulation:

A1 = 2.5*10^(-3);
A2 = A1;
Be = 1.0*10^9;
M = 15;
V1 = 1.0*10^(-3);
V2 = V1;
Cq = 0.67;
a = 0.05;
Cf = 10;
rho = 900;
x0 = 0.003;
Ps0 = 50*10^6;
P10 = 1*10^6;
P20 = 1*10^6;

p1 = Cq*a*x0/(2*sqrt(rho*(Ps0 -P10)))
q1 = Cq*a*sqrt(2/rho*(Ps0 - P10))

p2 = Cq*a*x0/(2*sqrt(rho*(P20)))
q2 = Cq*a*sqrt(2/rho*(P20))

Cl = 0.1; % Leak

A = [-Be*(p1+Cl)/V1 Be*Cl/V1 0 -Be*A1/V1; 
     -Be*Cl/V2 -Be*(p2-Cl)/V2 0 Be*A2/V2; 
     0 0 0 1;
     A1/M -A2/M 0 -Cf/M]

B = [Be*q1/V1 0; 
     Be*q2/V1 0;
     0 0;
     0 -1/M]

eig(A) % Unstable pole     

sys = ss(0, A, B);
step(sys);

y1 is $P_1$, y2 is $P_2$, y3 is $y$ and y4 is $\dot{y}$.

So if we add a stiffness $K = 5$ in to our model:

$$\begin{bmatrix} \dot{P_1}\\ \dot{P_2}\\ \dot{y_1}\\ \dot{y_2} \end{bmatrix} = \begin{bmatrix} -\frac{\beta (p_1 + C_l)}{V_1} & \frac{\beta C_l}{V_1} & 0 & -\frac{\beta A_1}{V_1}\\ \frac{\beta C_l}{V_2}& -\frac{\beta (p_2+C_l)}{V_2} &0 & \frac{\beta A_2}{V_2}\\ 0 & 0 &0 &1 \\ \frac{A_1}{M}& -\frac{A_2}{M} &-\frac{K}{M} & -\frac{C_f}{M} \end{bmatrix}\begin{bmatrix} P_1\\ P_2\\ y_1\\ y_2 \end{bmatrix} + \begin{bmatrix} \frac{\beta q_1}{V_1} & 0\\ -\frac{\beta q_2}{V_2} & 0\\ 0 & 0\\ 0 & -\frac{1}{M} \end{bmatrix}\begin{bmatrix} x\\ F_L \end{bmatrix} \\ \begin{bmatrix} y \end{bmatrix} = \begin{bmatrix} 0 &0 &1 &0 \end{bmatrix}\begin{bmatrix} P_1\\ P_2\\ y_1\\ y_2 \end{bmatrix}+\begin{bmatrix} 0 &0 \end{bmatrix}\begin{bmatrix} x\\ F_L \end{bmatrix}$$

We would have y3 to act like there is a limit of position:

Yes it is!

How did I found the solution to this? Answer: Equation (10.27) from the book Principles of Hydraulic Systems Design Second Edition.