Why is my state space model of a hydraulic servo system unstable?

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}$$

• engineering.stackexchange.com may be a better fit. Commented Mar 16, 2018 at 9:19
• @RodrigodeAzevedo But it's still math about state space. Commented Mar 16, 2018 at 9:50
• There are control engineers there, though: engineering.stackexchange.com/questions/tagged/… Commented Mar 16, 2018 at 10:01
• @RodrigodeAzevedo The question is on-topic here as well, so when the OP seems to object I won't migrate. If somebody who is very active at Engineering.SE also suggests migraion, I will reconsider. Blind migrations are not recommended. If I migrate I want to be sure that the receiving site approves. Having said all that I'm leaving the migration suggestion open. I simply don't know whether that would be a clearly better site for this. When in doubt, the OP's opinion gets extra weight. Commented Mar 16, 2018 at 11:11
• I think the unstability derives from the linearization itself. It may be eventually unstable since the linearization is valid only for a limited range of values. Is the system stable in a reasonable range? Commented Mar 16, 2018 at 12:39

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.

• Oh wow. Thank you for this post, Daniel. This is so helpful. Commented May 21, 2021 at 20:12
• :) :) :) :) :) :) Commented May 21, 2021 at 20:14