# 84. Multivariable pairing (RGA)¶

For a 2$:nbsphinx-math:times$2 system, we have 2 choices of pairing variables for distributed control:

 Diagonal Off-diagonal $\begin{split}G_{cd} = \left[\begin{array}{cc} G_{c1} & 0 \\ 0 & G_{c2} \end{array}\right]\end{split}$ $\begin{split}G_{co} = \left[\begin{array}{cc} 0 & G_{c2} \\ G_{c1} & 0 \end{array}\right]\end{split}$

Bristol developed the Relative Gain Array to determine good pairings based on only the plant transfer function matrix $$G_p$$. The elements of the RGA are defined as

$\lambda_{ij} \triangleq \frac{(\partial y_i/\partial u_j)_u}{(\partial y_i/\partial u_j)_y}= \frac{\text{open loop gain}}{\text{closed loop gain}}$

We could build $$\Lambda$$ by direct evaluation of the above derivatives near some point given a time-domain model, but if we already have a transfer function model, we can evaluate the steady-state gain matrix $$K$$ by using the final value theorem.

[1]:

import sympy
sympy.init_printing()

[2]:

s = sympy.Symbol('s')

[3]:

def fopdt(k, theta, tau):
return k*sympy.exp(-theta*s)/(tau*s + 1)


Using the system from example 16.5

[4]:

G_p = sympy.Matrix([[fopdt(-2, 1, 10), fopdt(1.5, 1, 1)],
[fopdt(1.5, 1, 1), fopdt(2, 1, 10)]])
G_p

[4]:

$$\left[\begin{matrix}- \frac{2 e^{- s}}{10 s + 1} & \frac{1.5 e^{- s}}{s + 1}\\\frac{1.5 e^{- s}}{s + 1} & \frac{2 e^{- s}}{10 s + 1}\end{matrix}\right]$$

Unfortunately sympy cannot calculate limits on matrix expressions

[5]:

#K = sympy.limit(G_p, s, 0)


But we can apply a function to the elements:

[6]:

def gain(G):
return sympy.limit(G, s, 0)

[7]:

K = G_p.applyfunc(gain)

[8]:

K

[8]:

$$\left[\begin{matrix}-2 & 1.5\\1.5 & 2\end{matrix}\right]$$

We can then calculate $$\Lambda = K \otimes H$$ where $$H = (K^{-1})^{T}$$:

[9]:

Lambda = K.multiply_elementwise(K.inv().transpose())
Lambda

[9]:

$$\left[\begin{matrix}0.64 & 0.36\\0.36 & 0.64\end{matrix}\right]$$

We can do the same calculation (faster) using numpy:

[10]:

import numpy

[11]:

def fopdt(k, theta, tau):
return k*numpy.exp(-theta*s)/(tau*s + 1)

[12]:

s = 0

[13]:

K = numpy.matrix([[fopdt(-2, 1, 10), fopdt(1.5, 1, 1)],
[fopdt(1.5, 1, 1), fopdt(2, 1, 10)]])


The .A attribute in matrices is the matrix as a numpy.array, which multiplies elementwise by default.

[14]:

K.A*K.I.T.A

[14]:

array([[0.64, 0.36],
[0.36, 0.64]])


The numpy developers recommended that you should use numpy.array instead of numpy.matrix as much as possible. I find this makes the notation harder to read:

[15]:

K = numpy.array([[fopdt(-2, 1, 10), fopdt(1.5, 1, 1)],
[fopdt(1.5, 1, 1), fopdt(2, 1, 10)]])

[16]:

K*numpy.linalg.inv(K).T

[16]:

array([[0.64, 0.36],
[0.36, 0.64]])


## 84.1. Simulation results¶

Let’s simulate this system to get an idea of how the control works out

[17]:

import tbcontrol
tbcontrol.expectversion('0.1.4')
from tbcontrol import blocksim

[18]:

import numpy

[19]:

N = 2

[20]:

G = {}

[21]:

gains = [[-2, 1.5],
[1.5, 2]]
taus = [[10, 1],
[1, 10]]
delays = [[1, 1],
[1, 1]]

[22]:

for inp in range(N):
for outp in range(N):
G[(outp, inp)] = blocksim.LTI(f"G_{outp}_{inp}", f"u_{inp}", f"yp_{inp}_{outp}",
gains[outp][inp], [1, taus[outp][inp]], delays[outp][inp])

[23]:

inputs = {'ysp_0': blocksim.step(),
'ysp_1': blocksim.step(starttime=50)}

[24]:

sums = {f'y_{outp}': [f"+yp_{inp}_{outp}" for inp in range(N)] for outp in range(N)}
for i in range(N):
sums[f'e_{i}'] = [f'+ysp_{i}', f'-y_{i}']

[25]:

sums

[25]:

{'y_0': ['+yp_0_0', '+yp_1_0'],
'y_1': ['+yp_0_1', '+yp_1_1'],
'e_0': ['+ysp_0', '-y_0'],
'e_1': ['+ysp_1', '-y_1']}

[26]:

import matplotlib.pyplot as plt
%matplotlib inline

[27]:

def simulate(auto1=True, K1=-1, tauI1=10, auto2=True, K2=0.5, tauI2=10):
controllers = {'Gc_0': blocksim.PI('Gc_0', 'e_0', 'u_0', K1, tauI1),
'Gc_1': blocksim.PI('Gc_1', 'e_1', 'u_1', K2, tauI2)}

controllers['Gc_0'].automatic = auto1
controllers['Gc_1'].automatic = auto2

ts = numpy.arange(0, 100, 0.125)

diagram = blocksim.Diagram(list(G.values()) + list(controllers.values()), sums, inputs)

result = diagram.simulate(ts)

#     plt.figure()
#     plt.plot(ts, result['u_0'])
#     plt.plot(ts, result['u_1'])

plt.figure()
plt.plot(ts, result['y_0'])
plt.plot(ts, result['ysp_0'])
plt.plot(ts, result['y_1'])
plt.plot(ts, result['ysp_1'])

[28]:

from ipywidgets import interact

[29]:

interact(simulate,
auto1=[True, False], K1=(-2., 0), tauI1=(1., 50),
auto2=[True, False], K2=(0., 2), tauI2=(1., 50))

[29]:

<function __main__.simulate(auto1=True, K1=-1, tauI1=10, auto2=True, K2=0.5, tauI2=10)>