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.
This conclusion is correct since all poles have a negative real part. However if the transfer function would have had more then one pole on the imaginary axis it would not have been possible to conclude whether the system is Lyapunov stable or not.
A more general approach would be to calculate the minimal state-space model. If its $A$ matrix has only eigenvalues with negative real then the system is Lyapunov, BIBO and asymptotically stable. Since it is a minimal realization the system would not be BIBO or asymptotically stable whenever the $A$ matrix has an eigenvalue along the imaginary axis, but the system still could be Lyapunov stable. Namely the system is still Lyapunov stable whenever the geometric multiplicity of the eigenvalues on the imaginary axis is at most one.
Best Answer
I am guessing that you are looking for the transfer function from $u$ to $y$, this would be consistent with current nomenclature.
Taking Laplace transforms gives $$ (s^2+2s) \hat{y_1} + s\hat{y_2} + \hat{u_1} = 0\\ (s-1)\hat{y_2} + \hat{u_2}-s \hat{u_1} = 0 $$ Solving algebraically gives $$\hat{y_1} = \frac{1-s-s^2}{s(s+2)(s-1)} \hat{u_1} + \frac{1}{s(s+2)(s-1)}\hat{u_2} \\ \hat{y_2} = \frac{s}{s-1} \hat{u_1} -\frac{1}{s-1} \hat{u_2} $$ from which all four transfer functions can be read off.