4.7. Observer Synthesis using Linear Matrix Inequalities#

Note: This is an advanced topic

4.7.1. References#

  1. Charkabarty, Ankush. “Guest Lecture: Exploiting Linear Matrix Inequalities In Control Systems Design”. (2015). Retrieved 29 March 2021.

  2. Boyd, Stephen, et al. Linear matrix inequalities in system and control theory. Society for industrial and applied mathematics, 1994.

  3. Caverly, Ryan James, and James Richard Forbes. “Lmi properties and applications in systems, stability, and control theory.” arXiv preprint arXiv:1903.08599 (2019).

4.7.2. Model Parameters#

4.7.3. Lyapunov Design#

Assuming no modeling error and ignoring disturbance inputs, the observer dynamics are described by

dedt=(ALC)e

where e=x^x is the difference between the estimated and process states. Given a symmetric positive definite matrix P, define the Lyapunov frunction V(e) as V(e)=eTPe.

dVdt=e˙TPe+eTPe˙=eT(ALC)TPe+eTP(ALC)e=eT(ATP+PACTLTPPLC)e

A sufficient condition for the global asympototic stability of observer is the left-hand side of this equation be negative for all e0. This will be true if and only if ATP+PACTLTPPLC is negative definite, i.e;,

eT(ATP+PACTLTPPLC)e<0e0

To provide some margin for robustness relative to model error, we will specify

dVdtγV

for some γ>0. When recast as a linear matrix inequality, we obtain

ATP+PACTLTPPLCγP

where the notation Q0 implies Q is negative definite.

As stated, given parameters A, C and γ, the task is to find a symmetric positive definite P0 and a matrix of observer gains L which satistify the condition above. As stated, this is a bilinear relationship due to the terms PL and CTLT appearing in the expression. A standard ‘trick’ is to introduce a new decision variable Y=PL to yield the linear matrix inequality

ATP+PACTYTYC+γP0

After finding a satisfactory solution P0 and Y, the desired observer gains are given by

L=P1Y

The next challenge is to perform the necessary calculations.

The first challenge is that the above relationship is homogeneous. In other words, for any scale factor α>0, if P and Y satisfy the relationship then so do αP and αY resulting in the same solution for L. Without loss of generality, we can specify a specific solution by adding the constraints

PI

4.7.4. CVXPY#

!pip install cvxpy
# parameter estimates.
alpha = 0.00016       # watts / (units P * percent U1)
P1 = 200              # P units
P2 = 100              # P units
CpH = 4.46            # heat capacity of the heater (J/deg C)
CpS = 0.819           # heat capacity of the sensor (J/deg C)
Ua = 0.050            # heat transfer coefficient from heater to environment
Ub = 0.021            # heat transfer coefficient from heater to sensor
Uc = 0.0335           # heat transfer coefficient between heaters
Tamb = 21             # ambient room temperature

# state space model
A = np.array([[-(Ua + Ub + Uc)/CpH, Ub/CpH, Uc/CpH, 0], 
              [Ub/CpS, -Ub/CpS, 0, 0],
              [Uc/CpH, 0, -(Ua + Ub + Uc)/CpH, Ub/CpH],
              [0, 0, Ub/CpS, -Ub/CpS]])

Bu = np.array([[alpha*P1/CpH, 0], [0, 0], [0, alpha*P2/CpH], [0, 0]])

Bd = np.array([[Ua/CpH], [0], [Ua/CpH], [0]])

C = np.array([[0, 1, 0, 0], [0, 0, 0, 1]])

n = A.shape[0]
p = C.shape[0]
import numpy as np
import cvxpy as cp

P = cp.Variable((n, n), PSD=True)
Y = cp.Variable((n, p))

gamma = 1/10
constraints = [P >> np.eye(n)]
constraints += [A.T@P + P@A - C.T@Y.T - Y@C + gamma*P << 0]

prob = cp.Problem(cp.Minimize(0), constraints)
prob.solve()
L = np.linalg.inv(P.value)@Y.value
print(L)
[[0.56726483 0.34537519]
 [0.22509201 0.09832281]
 [0.34537519 0.56726483]
 [0.09832281 0.22509201]]

4.7.5. H2 Optimal State Estimation#

Let’s now consider performance metrics for the observer. In particular, we assume our system is of the form

dxdt=Ax+Bddy=Cyxz=Czx

where d are unmeasured disturbances, y are process measurements, and z are process variables we are attempting to estimate. Constructing an estimator

dx^dt=Ax^L(y^y)+Bdd^y^=Cyx^

and defining the state error in the usual way as ex=x^x yields error dynamics given by

dedt=(ALCy)e+Bd(d^d)ez=Cze

where ez=Cze denotes the estimator error of interest. The design objective is to find values for L that minimize the impact of changes in disturbance d^d on ez.

[PA+ATPYCyCyTYTPBdBdTPI]0[PCzTCzZ]0Tr(Z)<ν
import numpy as np
import cvxpy as cp

P = cp.Variable((n, n), PSD=True)
Z = cp.Variable((p, p), PSD=True)
Y = cp.Variable((n, p))
nu = cp.Variable()

Cz = np.array([[1, 0, 0, 0], [0, 0, 1, 0]])

constraints = [cp.bmat([[A.T@P + P@A - Y@C - C.T@Y.T + np.eye(n), P@Bd], 
                       [Bd.T@P, -np.ones((1,1))]]) << 0]
constraints += [cp.bmat([[P, Cz.T], [Cz, Z]]) >> 0]
constraints += [cp.trace(Z) <= nu]
constraints += [nu >= 0]

prob = cp.Problem(cp.Minimize(nu), constraints)
prob.solve()
print(nu.value)
L = np.linalg.inv(P.value) @ Y.value
print(L)
0.01510843285193801
[[ 0.02450293 -0.00062794]
 [ 0.00531293  0.00279802]
 [-0.00062794  0.02450293]
 [ 0.00279802  0.00531293]]
import numpy as np
import cvxpy as cp

P = cp.Variable((5, 5), PSD=True)
Z = cp.Variable((1, 1), PSD=True)
Y = cp.Variable((5, 2))
nu = cp.Variable()


A_aug = np.vstack([np.hstack([A, Bd]), np.zeros([1, 5])])
Bu_aug = np.vstack([Bu, [[0, 0]]])
Bd_aug = np.vstack([np.zeros([4, 1]), [[1]]])
C_aug = np.hstack([C, np.zeros([2, 1])])

Cz_aug = np.array([[0, 0, 0, 0, 1]])


P = cp.Variable((5, 5), PSD=True)
Y = cp.Variable((5, 2))

gamma = 1/20
constraints = [P >> np.eye(5)]
constraints += [A_aug.T@P + P@A_aug - C_aug.T@Y.T - Y@C_aug + gamma*P << 0]

prob = cp.Problem(cp.Minimize(0), constraints)
prob.solve()
L = np.linalg.inv(P.value)@Y.value
print(L)

constraints = [cp.bmat([[A_aug.T@P + P@A_aug - Y@C_aug - C_aug.T@Y.T + np.eye(5), P@Bd_aug], 
                       [Bd_aug.T@P, -np.ones((1,1))]]) << 0]
constraints += [cp.bmat([[P, Cz_aug.T], [Cz_aug, Z]]) >> 0]
constraints += [cp.trace(Z) <= nu]
constraints += [nu >= 0]

prob = cp.Problem(cp.Minimize(0), constraints)
prob.solve()
print(P.value)
print(nu.value)
L = np.linalg.inv(P.value) @ Y.value
print(L)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-baf7d3ce8340> in <module>
      8 
      9 
---> 10 A_aug = np.vstack([np.hstack([A, Bd]), np.zeros([1, 5])])
     11 Bu_aug = np.vstack([Bu, [[0, 0]]])
     12 Bd_aug = np.vstack([np.zeros([4, 1]), [[1]]])

NameError: name 'A' is not defined

4.7.6. H Optimal State Estimation#

import numpy as np
import cvxpy as cp

P = cp.Variable((n, n), PSD=True)
Y = cp.Variable((4, 2))

constraint = cp.bmat([[A.T@P + P@A - Y@C - C.T@Y.T + np.eye(n), P@Bd], 
                       [Bd.T@P, -2*np.ones((1,1))]]) << 0
prob = cp.Problem(cp.Minimize(0), [constraint])
prob.solve()
L = np.linalg.inv(P.value) @ Y.value
print(L)
[[ 0.13413964  0.02750461]
 [ 0.21936083 -0.01039676]
 [ 0.02750461  0.13413964]
 [-0.01039676  0.21936083]]
ev, _ = np.linalg.eig(A - L@C)
-1/ev
array([21.53881154+20.40481958j, 21.53881154-20.40481958j,
       24.71082258+21.80171218j, 24.71082258-21.80171218j])
(repr(L))
'array([[ 0.0247742 , -0.00252384],\n       [-0.00184805,  0.00579824],\n       [-0.00252384,  0.0247742 ],\n       [ 0.00579824, -0.00184805]])'