Consider a hydropower plant with a dam. We want to control the flow through the dam gates in order to keep the amount of water at a desired level.
With $F_b(z) = \frac{S(z)}{R(z)}$ and $F_f(z) = \frac{T(z)}{R(z)}$, and using Mason's rule, we get that the closed-loop pulse-transfer function from command signal $u_c$ to output $y$ becomes $$G_c(z) = \frac{\frac{T(z)}{R(z)}\frac{b_u}{z-1}}{1 + \frac{S(z)}{R(z)} \frac{b_u}{(z-1)z}} = \frac{b_uzT(z)}{z(z-1)R(z) + b_uS(z)}.$$ The closed-loop transfer function from disturbance to output becomes $$G_{cv}(z) = \frac{\frac{b_v}{z-1}}{1 + \frac{S(z)}{R(z)} \frac{b_u}{(z-1)z}} = \frac{b_vzR(z)}{z(z-1)R(z) + b_uS(z)}.$$
The diophantine equation becomes $$z(z-1)R(z) + b_uS(z) = A_c(z)A_o(z)$$ We want to find the smallest order controller that can satisfy the Diophantine equation. Since the feedback controller is given by $$ F_b(z) = \frac{s_0z^n + s_1z^{n-1} + \cdots + s_n}{z^n + r_1z^{n-1} + \cdots + r_n}$$ and has $2\deg R + 1$ unknown parameters, and since we should choose the order of the Diphantine equation to be the same as the number of unknown parameters, we get $$ \deg \big((z(z-1)R(z) + b_uS(z)\big) = \deg R + 2 = 2\deg R + 1 \quad \Rightarrow \quad \deg R = n = 1.$$
The Diophantine equation thus becomes $$ z(z-1)(z+r_1) + b_u(s_0z+s_1) = (z-0.7)^2(z-0.5), $$ where $A_o(z) = z-0.5$ is the observer polynomial. Working out the expressions on both sides gives $$ z^3-(1-r_1)z^2 -r_1 z + b_us_0z + b_us_1 = (z^2 - 1.4z + 0.49)(z-0.5)$$ $$ z^3 -(1-r_1)z^2 +(b_us_0-r_1)z + b_us_1 = z^3 - (1.4+0.5)z^2 + (0.49+0.7)z -0.245$$ From the Diophantine equation we get the following equations in the unknowns \begin{align} z^2: &\quad 1-r_1 = 1.9\\ z^1: &\quad b_us_0 - r_1 = 1.19\\ z^0: &\quad b_us_1 = -0.245 \end{align} This is a linear system of equations in the unknown, and can be solved in many different ways. Here we see that with simple substitution we find \begin{align} r_1 &= 1-1.9 = -0.9\\ s_0 &= \frac{1}{b_u}(1.19+r_1) = \frac{0.29}{b_u}\\ s_1 &= -\frac{0.245}{b_u} \end{align}
We set $T(z) = t_0A_o(z)$ which gives the closed-loop pulse-transfer function $$G_c(z) = \frac{b_uzT(z)}{z(z-1)R(z) + b_uS(z)}= \frac{b_ut_0zA_o(z)}{A_c(z)A_o(z)} = \frac{b_u t_0z}{A_c(z)}$$ In order for this pulse-transfer function to have unit DC-gain (static gain) we must have $G_c(1) = 1$, or $$ \frac{b_ut_0}{A_c(1)} = 1. $$ The solution is $$ t_0 = \frac{A_c(1)}{b_u} = \frac{(1-0.7)^2}{b_u} = \frac{0.3^2}{b_u}. $$
In [1]:
import numpy as np
import sympy as sy
In [2]:
z = sy.symbols('z', real=False)
bu,r1,s0,s1 = sy.symbols('bu,r1,s0,s1', real=True)
pc,po = sy.symbols('pc,po', real=True) # Closed-loop pole and observer pole
In [17]:
# The polynomials
Ap = sy.Poly(z*(z-1), z)
Bp = sy.Poly(bu,z)
Rp = sy.Poly(z+r1, z)
Sp = sy.Poly(s0*z+s1, z)
Ac = sy.Poly((z-pc)**2, z)
Ao = sy.Poly(z-po, z)
In [18]:
# The diophantine eqn
dioph = Ap*Rp + Bp*Sp - Ac*Ao
# Form system of eqs from coefficients, then solve
dioph_coeffs = dioph.all_coeffs()
# Solve for r1, s0 and s1,
sol = sy.solve(dioph_coeffs, (r1,s0,s1))
print('r_1 = %s' % sol[r1])
print('s_0 = %s' % sol[s0])
print('s_1 = %s' % sol[s1])
In [19]:
# Substitute values for the desired closed-loop pole and observer pole
substitutions = [(pc, 0.7), (po, 0.5)]
print('r_1 = %s' % sol[r1].subs(substitutions))
print('s_0 = %s' % sol[s0].subs(substitutions))
print('s_1 = %s' % sol[s1].subs(substitutions))
In [20]:
# The forward controller
t0 = (Ac.eval(1)/Bp.eval(1))
print('t_0 = %s' % t0)
print('t_0 = %s' % t0.subs(substitutions))
Notice the solution for the controller denominator $$ R(z) = z+r_1 = z -2p_c -p_o + 1, $$ where $0\ge p_c<1$ is the desired closed-loop pole and $0 \ge p_o<1$ is the observer pole. Sketch in the $(p_c, p_o)$-plane the region which will give a stable controller $F_b(z) = \frac{S(z)}{R(z)}$!
In [30]:
import control
import control.matlab as cm
import matplotlib.pyplot as plt
In [33]:
sbs = [(bu, 1), (pc, 2.0/3.0), (po, 2.0/3.0)]
Rcoeffs = [1, float(sol[r1].subs(sbs))]
Scoeffs = [float(sol[s0].subs(sbs)), float(sol[s1].subs(sbs))]
Tcoeffs = float(t0.subs(sbs))*np.array([1, float(pc.subs(sbs))])
Acoeffs = [1, -1]
H = cm.tf(float(bu.subs(sbs)), Acoeffs, 1)
Ff = cm.tf(Tcoeffs, Rcoeffs, 1)
Fb = cm.tf(Scoeffs, Rcoeffs, 1)
Haa = cm.tf(1, [1, 0], 1) # The delay due to the anti-aliasing filter
Gc = cm.minreal(Ff*cm.feedback(H, Haa*Fb ))
Gcv = cm.feedback(H, Haa*Fb)
# Pulse trf fcn from command signal to control signal
Gcu = Ff*cm.feedback(1, H*Haa*Fb)
In [38]:
cm.pzmap(Fb)
Out[38]:
In [32]:
tvec = np.arange(40)
(t1, y1) = control.step_response(Gc,tvec)
plt.figure(figsize=(14,4))
plt.step(t1, y1[0])
plt.xlabel('k')
plt.ylabel('y')
plt.title('Output')
Out[32]:
In [40]:
(t1, y1) = control.step_response(Gcv,tvec)
plt.figure(figsize=(14,4))
plt.step(t1, y1[0])
plt.xlabel('k')
plt.ylabel('y')
plt.title('Output')
Out[40]:
In [ ]: