# Control of a hydropower dam
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. 
<p></p>
The system is a typical integrator, and is given by the difference equation
$$ y(k+1) = y(k) + b_uu(k) - b_vv(k), $$
where $x$ is the deviation of the water level from a reference level, $u$ is the change in the flow through the dam gates. A positive value of $u$ corresponds to less flow through the gates, relative to an operating point. The flow $v$ corresponds to changes in the flow in (from river) or out (through power plant). 
The pulse transfer function of the dam is thus $$H(z) = \frac{b_u}{z-1}.$$
We want to control the system using a two-degree-of-freedom controller, including an anti-aliasing filter modelled as a delay of one sampling period. This gives the block diagram <p></p>
The desired closed-loop system from the command signal $u_c$ to the output $y$ should have poles in $z=0.7$, and any observer poles should be chosen faster than the closed-loop poles, say in $z=0.5$.

## The closed-loop pulse-transfer functions
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
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}

## The feedforward 
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}. $$

## Verify by symbolic computer algebra


```python
import numpy as np
import sympy as sy
```


```python
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
```


```python
# 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)
```


```python
# 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])
```

    r_1 = -2*pc - po + 1
    s_0 = (pc**2 + 2*pc*po - 2*pc - po + 1)/bu
    s_1 = -pc**2*po/bu



```python
# 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))
```

    r_1 = -0.900000000000000
    s_0 = 0.29/bu
    s_1 = -0.245/bu



```python
# The forward controller
t0 = (Ac.eval(1)/Bp.eval(1))
print('t_0 = %s' % t0)
print('t_0 = %s' % t0.subs(substitutions))
```

    t_0 = (pc**2 - 2*pc + 1)/bu
    t_0 = 0.0900000000000001/bu


## Requirements on the closed-loop poles and observer poles in order to obtain stable controller
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)}$!

## Simulate a particular case
Let $b_u=1$, $p_c = p_o = \frac{2}{3}$. Analyze the closed-loop system by simulation


```python
import control
import control.matlab as cm
import matplotlib.pyplot as plt
```


```python
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)
```

    1 states have been removed from the model



```python
cm.pzmap(Fb)
```


```python
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')
```


```python
(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')
```


```python

```

cb8b1091e7b54ce77f4f78f5ff3bf19c9cced9b6
ipynb
Jupyter Notebook
polynomial-design/notebooks/RST-example-hydropower-dam.ipynb
kjartan-at-tec/mr2007-computerized-control
16e35f5007f53870eaf344eea1165507505ab4aa
MIT
2020-11-07T05:20:37.000Z
2020-12-22T09:46:13.000Z
polynomial-design/notebooks/RST-example-hydropower-dam.ipynb
alfkjartan/control-computarizado
5b9a3ae67602d131adf0b306f3ffce7a4914bf8e
MIT
2020-06-12T20:44:41.000Z
2020-06-12T20:49:00.000Z
polynomial-design/notebooks/RST-example-hydropower-dam.ipynb
alfkjartan/control-computarizado
5b9a3ae67602d131adf0b306f3ffce7a4914bf8e
MIT
2019-09-25T20:02:23.000Z
2019-09-25T20:02:23.000Z
Qwen/Qwen-72B
1. YES
2. YES