4.7. Observer Synthesis using Linear Matrix Inequalities#
Note: This is an advanced topic
4.7.1. References#
Charkabarty, Ankush. “Guest Lecture: Exploiting Linear Matrix Inequalities In Control Systems Design”. (2015). Retrieved 29 March 2021.
Boyd, Stephen, et al. Linear matrix inequalities in system and control theory. Society for industrial and applied mathematics, 1994.
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
where
A sufficient condition for the global asympototic stability of observer is the left-hand side of this equation be negative for all
To provide some margin for robustness relative to model error, we will specify
for some
where the notation
As stated, given parameters
After finding a satisfactory solution
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
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. Optimal State Estimation#
Let’s now consider performance metrics for the observer. In particular, we assume our system is of the form
where
and defining the state error in the usual way as
where
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. 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]])'