# CLOSED LOOP SYSTEMS

This notebook is a companion to the slides for understanding closed loop systems.

# Preliminaries


```python
!pip install -q control
!pip install -q tellurium
!pip install -q controlSBML

import control
import controlSBML as ctl
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sympy
import tellurium as te
```


```python
TIMES = ctl.makeSimulationTimes(0, 5, 500)
```

# Helpers


```python
def plotStep(tf, is_plot=True):
    """
    Plots a transfer function step response
    
    Parameters
    ----------
    tf: control.TransferFunction
    """
    result = control.forced_response(tf, T=TIMES, U=1)
    if is_plot:
        plt.plot(result.t, result.outputs)

# TESTS
Gs = control.TransferFunction([2], [1, 3])
plotStep(Gs, is_plot=False)
print("OK!")
```

    OK!


# Antimony Model

$\$S_1 \xrightarrow{k_1 S_1} S_2$,
$S_2 \xrightarrow{k_2 S_3} S_3$,
$S_3 \xrightarrow{k_3 S_3} S_2$,
$S_2 \xrightarrow{k_2 S_3} \emptyset$,


```python
# Constants
CONSTANT_DCT = {"k1": 1, "k2": 2, "k3": 3, "k4": 4}
s = sympy.Symbol("s")
REF = 10
```


```python
MODEL = """
$S1 -> S2; k1*$S1
J1: S2 -> S3; k2*S2
J2: S3 -> S2; k3*S3
J3: S2 -> ; k4*S2

k1 = 0
k2 = 0
k3 = 0
k4 = 0
$S1 = 10
S2 = 0
S3 = 0
S4 = 0
"""
RR = te.loada(MODEL)
for k in CONSTANT_DCT.keys():
    RR[k] = CONSTANT_DCT[k]
RR.simulate()
RR.plot()
```




We are going to control ``S3`` by changing ``S2``. That is, $U(s) = S_2  (s)$
and $Y(s) = S_3 (s)$.

# Simulating the Closed Loop System


```python
ctlsb = ctl.ControlSBML(RR,  input_names=["S2"], output_names=["S3"])
```


```python
sys = ctlsb.makeNonlinearIOSystem("sys")
```


```python
ctlsb.makeStateSpace()
```




\[
\left(\begin{array}{rll|rll}
-3\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&2\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\
\hline
1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\
\end{array}\right)
\]




```python
def controller_outfcn(_, __, u, param_dct):
    k = param_dct.get("k", 2)
    ref = param_dct.get("ref", REF)
    return k*([REF] - u[0])
#
controller = control.NonlinearIOSystem(
    None, controller_outfcn, inputs=['in'], outputs=['out'],
    name="controller")
```


```python
# Create the closed loop system
closed_loop = control.interconnect(
  [sys, controller],       # systems
  connections=[
    ['sys.S2', 'controller.out'],
    ['controller.in',  'sys.S3'],
  ],
  inplist=["controller.in"],
  outlist=["sys.S2", "sys.S3"],
)
```


```python
param_dct = {"k": 5, "ref": REF}
result = control.input_output_response(closed_loop, T=TIMES, params=param_dct)
plt.plot(TIMES, result.outputs[0], TIMES, result.outputs[1])
plt.legend(["S2", "S3"])
```


```python
result.outputs[1][-2:]
```




    array([7.68988725, 7.68994905])




```python
type(control.feedback(controller, sys))
```




    control.iosys.InterconnectedSystem



The output is biased! We wanted $S_3$ to be 10, but it's 7.6.

# Analyzing the Cause of Bias

Using direct calculation, I have calculated the transfer function $G(s) = \frac{k_2}{s + k_3}$.


```python
numerator = CONSTANT_DCT["k2"]
denominator = [1, CONSTANT_DCT["k3"]]
tf = control.TransferFunction(numerator, denominator)
tf
```




$$\frac{2}{s + 3}$$




```python
G = 2/(s + 3)
G
```




$\displaystyle \frac{2}{s + 3}$




```python
C = sympy.symbols("k")
```

We calculate the transfer function for the closed loop system, $H(s) = \frac{Y(s)}{R(s)}$.

# Constructing Transfer Functions



In the following, we use $G(s) = \frac{2}{s + 3}$ and $C(s) = k$ (proportional control).


```python
Gs = control.TransferFunction([2], [1, 3])
Gs
```




$$\frac{2}{s + 3}$$




```python
plotStep(Gs)
```

The original system, $G(s)$, settles around time 1.


```python
k = 5
Cs = control.TransferFunction([k], [1])
Cs.dcgain()
```




    5.0



For the filter, we use a low pass filter.

$F(s) = \frac{k_f}{s + k_f}$. Note that the DC Gain of $F(s)$ is 1.


```python
kf = 5
Fs = control.TransferFunction([kf], [1, kf])
Fs.dcgain()
```




    1.0



We use the above transfer functions to construct the various closed loop transfer functions.
Note that you can multiply and divide transfer function objects.


```python
HRYs = control.minreal(Cs*Gs/(1 + Cs*Gs))
HRYs
```

    1 states have been removed from the model





$$\frac{10}{s + 13}$$




```python
HRYs.pole()
```




    array([-13.])




```python
HRYs.dcgain()
```




    0.7692307692307693




```python
plotStep(HRYs)
```

The closed loop system settles around time 0.5.

Now we include the filter


```python
HRYs = control.minreal(Cs*Gs/(1 + Cs*Gs*Fs))
HRYs
```

    1 states have been removed from the model





$$\frac{10 s + 50}{s^2 + 8 s + 65}$$




```python
HRYs.dcgain()
```




    0.7692307692307698




```python
plotStep(HRYs)
```


```python
HREs = control.minreal(1/(1 + Cs*Gs*Fs))
HREs
```

    0 states have been removed from the model





$$\frac{s^2 + 8 s + 15}{s^2 + 8 s + 65}$$



We moved the pole of the original system from -3 to -13.

e464b598adaf8cd7801b72d2de4976739046d2df
ipynb
Jupyter Notebook
Lecture_13_14-Closed-Loop-Systems/Closed-Loop-Systems.ipynb
joseph-hellerstein/advanced-controls-lectures
dc43f6c3517616da3b0ea7c93192d911414ee202
MIT
Lecture_13_14-Closed-Loop-Systems/Closed-Loop-Systems.ipynb
joseph-hellerstein/advanced-controls-lectures
dc43f6c3517616da3b0ea7c93192d911414ee202
MIT
Lecture_13_14-Closed-Loop-Systems/Closed-Loop-Systems.ipynb
joseph-hellerstein/advanced-controls-lectures
dc43f6c3517616da3b0ea7c93192d911414ee202
MIT
Qwen/Qwen-72B
1. YES
2. YES