{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "*This notebook contains material from [cbe30338-2021](https://jckantor.github.io/cbe30338-2021);\n", "content is available [on Github](https://github.com/jckantor/cbe30338-2021.git).*\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "< [4.6 Lab Assignment 6: Anomaly Detection](https://jckantor.github.io/cbe30338-2021/04.06-Lab-Assignment-Anomaly-Detection.html) | [Contents](toc.html) | [Tag Index](tag_index.html) | [4.8 Application of Luenberger Observers to Environmental Modeling of Rivers](https://jckantor.github.io/cbe30338-2021/04.08-State-Estimation-Envrironmental-Application.html) >

\"Open

\"Download\"" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 1, "link": "[4.7 Observer Synthesis using Linear Matrix Inequalities](https://jckantor.github.io/cbe30338-2021/04.07-Observer-Synthesis-LMI.html#4.7-Observer-Synthesis-using-Linear-Matrix-Inequalities)", "section": "4.7 Observer Synthesis using Linear Matrix Inequalities" } }, "source": [ "# 4.7 Observer Synthesis using Linear Matrix Inequalities\n", "\n", "Note: This is an advanced topic" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[4.7.1 References](https://jckantor.github.io/cbe30338-2021/04.07-Observer-Synthesis-LMI.html#4.7.1-References)", "section": "4.7.1 References" } }, "source": [ "## 4.7.1 References\n", "\n", "1. Charkabarty, Ankush. [\"Guest Lecture: Exploiting Linear Matrix Inequalities In\n", "Control Systems Design\"](https://engineering.utsa.edu/ataha/wp-content/uploads/sites/38/2017/10/EE5243_Module10.pdf). (2015). Retrieved 29 March 2021.\n", "\n", "2. Boyd, Stephen, et al. [Linear matrix inequalities in system and control theory.](https://web.stanford.edu/~boyd/lmibook/index.html) Society for industrial and applied mathematics, 1994.\n", "\n", "3. Caverly, Ryan James, and James Richard Forbes. [\"Lmi properties and applications in systems, stability, and control theory.\"](https://arxiv.org/abs/1903.08599) arXiv preprint arXiv:1903.08599 (2019). \n" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[4.7.2 Model Parameters](https://jckantor.github.io/cbe30338-2021/04.07-Observer-Synthesis-LMI.html#4.7.2-Model-Parameters)", "section": "4.7.2 Model Parameters" } }, "source": [ "## 4.7.2 Model Parameters" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[4.7.3 Lyapunov Design](https://jckantor.github.io/cbe30338-2021/04.07-Observer-Synthesis-LMI.html#4.7.3-Lyapunov-Design)", "section": "4.7.3 Lyapunov Design" } }, "source": [ "## 4.7.3 Lyapunov Design\n", "\n", "Assuming no modeling error and ignoring disturbance inputs, the observer dynamics are described by\n", "\n", "$$\\frac{de}{dt} = (A - LC)e$$\n", "\n", "where $e = \\hat{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) = e^TPe$. \n", "\n", "\\begin{align}\n", "\\frac{dV}{dt} & = \\dot{e}^TP e + e^TP\\dot{e} \\\\\n", "& = e^T(A - LC)^T Pe + e^T P (A - LC) e \\\\ \n", "& = e^T(A^TP + PA - C^TL^TP - P L C)e\\\\\n", "\\end{align}\n", "\n", "A sufficient condition for the global asympototic stability of observer is the left-hand side of this equation be negative for all $e \\ne 0$. This will be true if and only if $A^TP + PA - C^TL^TP - P L C$ is negative definite, i.e;,\n", "\n", "$$e^T(A^TP + PA - C^TL^TP - P L C)e < 0 \\qquad\\forall e \\ne 0$$\n", "\n", "To provide some margin for robustness relative to model error, we will specify\n", "\n", "\\begin{align}\n", "\\frac{dV}{dt} \\leq -\\gamma V\n", "\\end{align}\n", "\n", "for some $\\gamma > 0$. When recast as a linear matrix inequality, we obtain\n", "\n", "$$A^TP + PA - C^TL^TP - P L C \\prec -\\gamma P$$\n", "\n", "where the notation $Q \\prec 0$ implies $Q$ is negative definite.\n", "\n", "As stated, given parameters $A$, $C$ and $\\gamma$, the task is to find a symmetric positive definite $P \\succ 0$ 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 $C^TL^T$ appearing in the expression. A standard 'trick' is to introduce a new decision variable $Y = PL$ to yield the linear matrix inequality\n", "\n", "$$A^TP + PA - C^TY^T - Y C + \\gamma P \\prec 0$$\n", "\n", "After finding a satisfactory solution $P \\succ 0$ and $Y$, the desired observer gains are given by\n", "\n", "$$L = P^{-1} Y$$\n", "\n", "The next challenge is to perform the necessary calculations. \n", "\n", "The first challenge is that the above relationship is homogeneous. In other words, for any scale factor $\\alpha > 0$, if $P$ and $Y$ satisfy the relationship then so do $\\alpha P$ and $\\alpha Y$ resulting in the same solution for $L$. Without loss of generality, we can specify a specific solution by adding the constraints\n", "\n", "$$P \\succ I$$" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[4.7.4 CVXPY](https://jckantor.github.io/cbe30338-2021/04.07-Observer-Synthesis-LMI.html#4.7.4-CVXPY)", "section": "4.7.4 CVXPY" } }, "source": [ "## 4.7.4 CVXPY\n", "\n", " !pip install cvxpy" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "nbpages": { "level": 2, "link": "[4.7.4 CVXPY](https://jckantor.github.io/cbe30338-2021/04.07-Observer-Synthesis-LMI.html#4.7.4-CVXPY)", "section": "4.7.4 CVXPY" } }, "outputs": [], "source": [ "# parameter estimates.\n", "alpha = 0.00016 # watts / (units P * percent U1)\n", "P1 = 200 # P units\n", "P2 = 100 # P units\n", "CpH = 4.46 # heat capacity of the heater (J/deg C)\n", "CpS = 0.819 # heat capacity of the sensor (J/deg C)\n", "Ua = 0.050 # heat transfer coefficient from heater to environment\n", "Ub = 0.021 # heat transfer coefficient from heater to sensor\n", "Uc = 0.0335 # heat transfer coefficient between heaters\n", "Tamb = 21 # ambient room temperature\n", "\n", "# state space model\n", "A = np.array([[-(Ua + Ub + Uc)/CpH, Ub/CpH, Uc/CpH, 0], \n", " [Ub/CpS, -Ub/CpS, 0, 0],\n", " [Uc/CpH, 0, -(Ua + Ub + Uc)/CpH, Ub/CpH],\n", " [0, 0, Ub/CpS, -Ub/CpS]])\n", "\n", "Bu = np.array([[alpha*P1/CpH, 0], [0, 0], [0, alpha*P2/CpH], [0, 0]])\n", "\n", "Bd = np.array([[Ua/CpH], [0], [Ua/CpH], [0]])\n", "\n", "C = np.array([[0, 1, 0, 0], [0, 0, 0, 1]])\n", "\n", "n = A.shape[0]\n", "p = C.shape[0]" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "nbpages": { "level": 2, "link": "[4.7.4 CVXPY](https://jckantor.github.io/cbe30338-2021/04.07-Observer-Synthesis-LMI.html#4.7.4-CVXPY)", "section": "4.7.4 CVXPY" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0.56726483 0.34537519]\n", " [0.22509201 0.09832281]\n", " [0.34537519 0.56726483]\n", " [0.09832281 0.22509201]]\n" ] } ], "source": [ "import numpy as np\n", "import cvxpy as cp\n", "\n", "P = cp.Variable((n, n), PSD=True)\n", "Y = cp.Variable((n, p))\n", "\n", "gamma = 1/10\n", "constraints = [P >> np.eye(n)]\n", "constraints += [A.T@P + P@A - C.T@Y.T - Y@C + gamma*P << 0]\n", "\n", "prob = cp.Problem(cp.Minimize(0), constraints)\n", "prob.solve()\n", "L = np.linalg.inv(P.value)@Y.value\n", "print(L)" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[4.7.5 $\\cal{H}_2$ Optimal State Estimation](https://jckantor.github.io/cbe30338-2021/04.07-Observer-Synthesis-LMI.html#4.7.5-$\\cal{H}_2$-Optimal-State-Estimation)", "section": "4.7.5 $\\cal{H}_2$ Optimal State Estimation" } }, "source": [ "## 4.7.5 $\\cal{H}_2$ Optimal State Estimation\n", "\n", "Let's now consider performance metrics for the observer. In particular, we assume our system is of the form\n", "\n", "\\begin{align}\n", "\\frac{dx}{dt} & = A x + B_d d \\\\\n", "y & = C_y x \\\\\n", "z & = C_z x\n", "\\end{align}\n", "\n", "where $d$ are unmeasured disturbances, $y$ are process measurements, and $z$ are process variables we are attempting to estimate. Constructing an estimator\n", "\n", "\\begin{align}\n", "\\frac{d\\hat{x}}{dt} & = A\\hat{x} - L(\\hat{y} - y) + B_d\\hat{d}\\\\\n", "\\hat{y} & = C_y\\hat{x}\n", "\\end{align}\n", "\n", "and defining the state error in the usual way as $e_x = \\hat{x} - x$ yields error dynamics given by\n", "\n", "\\begin{align}\n", "\\frac{de}{dt} & = (A - LC_y) e + B_d(\\hat{d} - d) \\\\\n", "e_z & = C_z e\n", "\\end{align}\n", "\n", "where $e_z = C_z e$ denotes the estimator error of interest. The design objective is to find values for $L$ that minimize the impact of changes in disturbance $\\hat{d} - d$ on $e_z$.\n", "\n", "\\begin{align}\n", "\\begin{bmatrix} PA + A^TP - YC_y - C_y^TY^T & PB_d \\\\ B_d^T P & -I\\end{bmatrix} & \\prec 0 \\\\\n", "\\begin{bmatrix} P & C_z^T \\\\ C_z & Z\\end{bmatrix} & \\succ 0 \\\\\n", "Tr(Z) & < \\nu\n", "\\end{align}" ] }, { "cell_type": "code", "execution_count": 268, "metadata": { "nbpages": { "level": 2, "link": "[4.7.5 $\\cal{H}_2$ Optimal State Estimation](https://jckantor.github.io/cbe30338-2021/04.07-Observer-Synthesis-LMI.html#4.7.5-$\\cal{H}_2$-Optimal-State-Estimation)", "section": "4.7.5 $\\cal{H}_2$ Optimal State Estimation" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.01510843285193801\n", "[[ 0.02450293 -0.00062794]\n", " [ 0.00531293 0.00279802]\n", " [-0.00062794 0.02450293]\n", " [ 0.00279802 0.00531293]]\n" ] } ], "source": [ "import numpy as np\n", "import cvxpy as cp\n", "\n", "P = cp.Variable((n, n), PSD=True)\n", "Z = cp.Variable((p, p), PSD=True)\n", "Y = cp.Variable((n, p))\n", "nu = cp.Variable()\n", "\n", "Cz = np.array([[1, 0, 0, 0], [0, 0, 1, 0]])\n", "\n", "constraints = [cp.bmat([[A.T@P + P@A - Y@C - C.T@Y.T + np.eye(n), P@Bd], \n", " [Bd.T@P, -np.ones((1,1))]]) << 0]\n", "constraints += [cp.bmat([[P, Cz.T], [Cz, Z]]) >> 0]\n", "constraints += [cp.trace(Z) <= nu]\n", "constraints += [nu >= 0]\n", "\n", "prob = cp.Problem(cp.Minimize(nu), constraints)\n", "prob.solve()\n", "print(nu.value)\n", "L = np.linalg.inv(P.value) @ Y.value\n", "print(L)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "nbpages": { "level": 2, "link": "[4.7.5 $\\cal{H}_2$ Optimal State Estimation](https://jckantor.github.io/cbe30338-2021/04.07-Observer-Synthesis-LMI.html#4.7.5-$\\cal{H}_2$-Optimal-State-Estimation)", "section": "4.7.5 $\\cal{H}_2$ Optimal State Estimation" } }, "outputs": [ { "ename": "NameError", "evalue": "name 'A' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 10\u001b[0;31m \u001b[0mA_aug\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvstack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhstack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mBd\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m5\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 11\u001b[0m \u001b[0mBu_aug\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvstack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mBu\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 12\u001b[0m \u001b[0mBd_aug\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvstack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'A' is not defined" ] } ], "source": [ "import numpy as np\n", "import cvxpy as cp\n", "\n", "P = cp.Variable((5, 5), PSD=True)\n", "Z = cp.Variable((1, 1), PSD=True)\n", "Y = cp.Variable((5, 2))\n", "nu = cp.Variable()\n", "\n", "\n", "A_aug = np.vstack([np.hstack([A, Bd]), np.zeros([1, 5])])\n", "Bu_aug = np.vstack([Bu, [[0, 0]]])\n", "Bd_aug = np.vstack([np.zeros([4, 1]), [[1]]])\n", "C_aug = np.hstack([C, np.zeros([2, 1])])\n", "\n", "Cz_aug = np.array([[0, 0, 0, 0, 1]])\n", "\n", "\n", "P = cp.Variable((5, 5), PSD=True)\n", "Y = cp.Variable((5, 2))\n", "\n", "gamma = 1/20\n", "constraints = [P >> np.eye(5)]\n", "constraints += [A_aug.T@P + P@A_aug - C_aug.T@Y.T - Y@C_aug + gamma*P << 0]\n", "\n", "prob = cp.Problem(cp.Minimize(0), constraints)\n", "prob.solve()\n", "L = np.linalg.inv(P.value)@Y.value\n", "print(L)\n", "\n", "constraints = [cp.bmat([[A_aug.T@P + P@A_aug - Y@C_aug - C_aug.T@Y.T + np.eye(5), P@Bd_aug], \n", " [Bd_aug.T@P, -np.ones((1,1))]]) << 0]\n", "constraints += [cp.bmat([[P, Cz_aug.T], [Cz_aug, Z]]) >> 0]\n", "constraints += [cp.trace(Z) <= nu]\n", "constraints += [nu >= 0]\n", "\n", "prob = cp.Problem(cp.Minimize(0), constraints)\n", "prob.solve()\n", "print(P.value)\n", "print(nu.value)\n", "L = np.linalg.inv(P.value) @ Y.value\n", "print(L)\n" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[4.7.6 $\\cal{H}_{\\infty}$ Optimal State Estimation](https://jckantor.github.io/cbe30338-2021/04.07-Observer-Synthesis-LMI.html#4.7.6-$\\cal{H}_{\\infty}$-Optimal-State-Estimation)", "section": "4.7.6 $\\cal{H}_{\\infty}$ Optimal State Estimation" } }, "source": [ "## 4.7.6 $\\cal{H}_{\\infty}$ Optimal State Estimation" ] }, { "cell_type": "code", "execution_count": 159, "metadata": { "nbpages": { "level": 2, "link": "[4.7.6 $\\cal{H}_{\\infty}$ Optimal State Estimation](https://jckantor.github.io/cbe30338-2021/04.07-Observer-Synthesis-LMI.html#4.7.6-$\\cal{H}_{\\infty}$-Optimal-State-Estimation)", "section": "4.7.6 $\\cal{H}_{\\infty}$ Optimal State Estimation" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 0.13413964 0.02750461]\n", " [ 0.21936083 -0.01039676]\n", " [ 0.02750461 0.13413964]\n", " [-0.01039676 0.21936083]]\n" ] } ], "source": [ "import numpy as np\n", "import cvxpy as cp\n", "\n", "P = cp.Variable((n, n), PSD=True)\n", "Y = cp.Variable((4, 2))\n", "\n", "constraint = cp.bmat([[A.T@P + P@A - Y@C - C.T@Y.T + np.eye(n), P@Bd], \n", " [Bd.T@P, -2*np.ones((1,1))]]) << 0\n", "prob = cp.Problem(cp.Minimize(0), [constraint])\n", "prob.solve()\n", "L = np.linalg.inv(P.value) @ Y.value\n", "print(L)" ] }, { "cell_type": "code", "execution_count": 191, "metadata": { "nbpages": { "level": 2, "link": "[4.7.6 $\\cal{H}_{\\infty}$ Optimal State Estimation](https://jckantor.github.io/cbe30338-2021/04.07-Observer-Synthesis-LMI.html#4.7.6-$\\cal{H}_{\\infty}$-Optimal-State-Estimation)", "section": "4.7.6 $\\cal{H}_{\\infty}$ Optimal State Estimation" } }, "outputs": [ { "data": { "text/plain": [ "array([21.53881154+20.40481958j, 21.53881154-20.40481958j,\n", " 24.71082258+21.80171218j, 24.71082258-21.80171218j])" ] }, "execution_count": 191, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ev, _ = np.linalg.eig(A - L@C)\n", "-1/ev" ] }, { "cell_type": "code", "execution_count": 194, "metadata": { "nbpages": { "level": 2, "link": "[4.7.6 $\\cal{H}_{\\infty}$ Optimal State Estimation](https://jckantor.github.io/cbe30338-2021/04.07-Observer-Synthesis-LMI.html#4.7.6-$\\cal{H}_{\\infty}$-Optimal-State-Estimation)", "section": "4.7.6 $\\cal{H}_{\\infty}$ Optimal State Estimation" } }, "outputs": [ { "data": { "text/plain": [ "'array([[ 0.0247742 , -0.00252384],\\n [-0.00184805, 0.00579824],\\n [-0.00252384, 0.0247742 ],\\n [ 0.00579824, -0.00184805]])'" ] }, "execution_count": 194, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(repr(L))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 2, "link": "[4.7.6 $\\cal{H}_{\\infty}$ Optimal State Estimation](https://jckantor.github.io/cbe30338-2021/04.07-Observer-Synthesis-LMI.html#4.7.6-$\\cal{H}_{\\infty}$-Optimal-State-Estimation)", "section": "4.7.6 $\\cal{H}_{\\infty}$ Optimal State Estimation" } }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "< [4.6 Lab Assignment 6: Anomaly Detection](https://jckantor.github.io/cbe30338-2021/04.06-Lab-Assignment-Anomaly-Detection.html) | [Contents](toc.html) | [Tag Index](tag_index.html) | [4.8 Application of Luenberger Observers to Environmental Modeling of Rivers](https://jckantor.github.io/cbe30338-2021/04.08-State-Estimation-Envrironmental-Application.html) >

\"Open

\"Download\"" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 4 }