{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 线性系统状态空间模型\n", "\n", "![关系](https://github.com/w407022008/All-of-Notes/blob/master/System_Science_and_Control_Engineering/state_space_model.png?raw=true)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## - 状态空间模型 State Space\n", "\n", "通常情况下，我们把一个线性系统描写作：\n", "$$\\dot X = AX + BU$$\n", "$$Y = CX + DU$$\n", "当中，矩阵**A**被称作**系统矩阵(System matrix)**，矩阵**B**被称作**输入矩阵(Input matrix)**，矩阵**C**被称作**输出矩阵(Output matrix)**，矩阵**D**被称作**前馈矩阵(Feedthrough matrix)**，向量**X**被称作**状态变量(State variable)**，向量**U**被称作**输入变量(Input)**，向量**Y**被称作**输出变量(Output)**。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## - 物理方程描述\n", "\n", "状态空间模型本身就是一个形似常微分方程的模型。在线性空间里，我们可以把状态变量**X**升维成为一个向量，如：$X=[X,\\dot X, \\ddot X]^T$，那么这个状态方程可以被写作：$$\\left ( \\begin{matrix} \\dot X \\\\ \\ddot X \\\\ \\dddot X \\end{matrix} \\right ) = \\left [\\begin{matrix} 0&1&0\\\\ 0&0&1\\\\ a_0&a_1&a_2 \\end{matrix}\\right] \\left(\\begin{matrix} X \\\\ \\dot X \\\\ \\ddot X \\end{matrix} \\right ) + \\left [\\begin{matrix} 0 \\\\ 0 \\\\ 1\\end{matrix}\\right]U$$ 于是，容易发现，他的物理模型是：$$\\dddot X = a_0X + a_1\\dot X + a_2\\ddot X + U$$\n", "\n", "在现实生活中，有很多的例子，他的物理方程是由某个高阶常微分方程构建的。对于这一类物理现象的描述，我们抽象出一个**数学模型：**\n", "$$y^{(n)} +a_1y^{(n-1)} + \\cdots +a_ny^{(0)}=b_0u^{(n)} +b_1u^{(n-1)} +\\cdots +b_nu^{(0)}$$\n", "我们令$$\\begin{cases} \\beta_0=b_0 \\\\ \\beta_1=b_1-a_1\\beta_0 \\\\ \\beta_2 = b_2-a_2\\beta_0-a_1\\beta_1 \\\\ \\vdots \\end{cases}且 \\begin{cases} x_1 =y -\\beta_0u \\\\ x_2 = \\dot y -\\beta_1u - \\beta_0\\dot u \\\\ x_3 = \\ddot y -\\beta_2u -\\beta_1 \\dot u -\\beta_0 \\ddot u \\\\ \\vdots \\\\ x_n = y^{(n-1)} -\\beta_{n-1}u -\\beta_{n-2} \\dot u - \\cdots -\\beta_0u^{(n-1)}\\end{cases}$$ 于是，有：$$\\begin{cases} \\dot X=\\begin{bmatrix} 0&1&0&\\cdots&0 \\\\ 0&0&1&\\cdots&0 \\\\ \\vdots&\\vdots&\\vdots&\\ddots&\\vdots \\\\ -a_n&-a_{n-1}&-a_{n-2}&\\cdots&-a_1\\end{bmatrix}X+\\begin{bmatrix} \\beta_1\\\\ \\beta_2\\\\ \\beta_3\\\\ \\vdots\\\\ \\beta_n \\end{bmatrix}U \\\\ Y=\\begin{bmatrix}1&0&0&\\cdots&0\\end{bmatrix}X+\\beta_0U \\end{cases}$$\n", "如果，**输入变量不含高阶导数项**，那么数学模型写作：$$y^{(n)} +a_1y^{(n-1)} + \\cdots +a_ny^{(0)}=bu$$ 状态空间模型写作：$$\\begin{pmatrix} \\dot y \\\\ \\ddot y \\\\ \\vdots \\\\ y^{(n)}\\end{pmatrix} = \\begin{bmatrix} 0&1&0&\\cdots&0 \\\\ 0&0&1&\\cdots&0 \\\\ \\vdots&\\vdots&\\vdots&\\ddots&\\vdots \\\\ -a_n&-a_{n-1}&-a_{n-2}&\\cdots&-a_1 \\end{bmatrix} \\begin{pmatrix} y \\\\ \\dot y\\\\ \\vdots \\\\ y^{(n-1)}\\end{pmatrix}+ \\begin{bmatrix} 0 \\\\ 0 \\\\ \\vdots \\\\ b\\end{bmatrix}u$$ \n", "\n", "**显然，不含输入导数项的物理方程更容易被用来构建状态空间模型**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example\n", "\n", "对于高阶微分方程：$\\dddot y + 3\\ddot y + 5\\dot y + 7y = 2\\dot u + 4u$，有：$$\\begin{cases}A=\\begin{bmatrix} 0&1&0\\\\0&0&1\\\\-7&-5&-3\\end{bmatrix} B=\\begin{bmatrix} 0\\\\2\\\\-2\\end{bmatrix}\\\\C=\\begin{bmatrix}1&0&0\\end{bmatrix} \\end{cases}$$" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "ans =\n", " \n", " 2 s + 4\n", " ---------------------\n", " s^3 + 3 s^2 + 5 s + 7\n", " \n", "Continuous-time transfer function.\n", "\n", "\n" ] } ], "source": [ "A=[0 1 0; 0 0 1; -7 -5 -3];\n", "B=[0;2;-2];\n", "C=[1 0 0];\n", "sys=ss(A,B,C,[]);% state space model\n", "tf(sys) % transfer function" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## - 输入输出传递函数模型\n", "\n", "在线性空间里，我们对物理方程的输入输出项取拉普拉斯变换，从而等到频域内输入输出传递函数：$$G(s)=\\frac{Y(s)}{U(s)}=\\frac{b_0s^n+b_1s^{n-1}+\\cdots+b_ns}{a_0s^n+a_1s^{n-1}+\\cdots+a_ns}=G_{\\Sigma(A,B,C)}(s)+\\frac{b_0}{a_0}$$\n", "于是，我们称方程$a_0s^n+a_1s^{n-1}+\\cdots+a_ns=0$为系统的**特征方程**，方程的根被称为**特征值**或**极点(Pole point)**，相应地，方程$b_0s^n+b_1s^{n-1}+\\cdots+b_ns=0$的根被称为**零点(Zero point)**。 \n", "当分母被写作$s^2+a_1s+a_2$时，有$\\omega=\\sqrt{a_2},\\xi=\\frac{a_1}{2\\omega}$，$\\omega$被称作固有频率，$\\xi$被称作阻尼比。$\\xi=1$表示临界阻尼，当$\\xi>1$时，系统有大阻尼比。 \n", "![](https://upload.wikimedia.org/wikipedia/commons/thumb/9/94/2nd_Order_Damping_Ratios.svg/800px-2nd_Order_Damping_Ratios.svg.png?1539979741329)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- 如果传递函数$G_{\\Sigma(A,B,C)}(s)$中，**极点互异（不重根）**时，传递函数可写作：\n", "$$G(s)=\\frac{k_1}{s-s_1}+\\frac{k_2}{s-s_2}+\\cdots+\\frac{k_n}{s-s_n}+D$$\n", "令：$\\begin{cases} X_i(s)=\\frac{1}{s-s_i}U(s) \\Rightarrow \\dot x_i=s_ix_i + u \\\\ Y(s)=\\sum k_iX_i(s)+Du \\Rightarrow y=k_1x_1+k_2x_2+\\cdots+k_nx_n+Du\\end{cases}$ \n", "$$\\Rightarrow \\begin{cases}\\dot X=\\begin{bmatrix}s_1&&&\\\\&s_2&&\\\\&&\\ddots&\\\\&&&s_n\\end{bmatrix}X+\\begin{bmatrix}1\\\\1\\\\ \\vdots\\\\1\\end{bmatrix}u \\\\ Y=\\begin{bmatrix}k_1&k_2&\\cdots&k_n\\end{bmatrix}X+Du\\end{cases}$$ \n", "\n", "- 如果**重极点**存在，如：$G(s)=\\frac{k_{11}}{(s-s_1)^3}+\\frac{k_{12}}{(s-s_1)^3}+\\frac{k_{13}}{(s-s_1)^3}+\\frac{k_{41}}{(s-s_4)^3}+\\frac{k_{51}}{(s-s_5)^3}+\\frac{k_{52}}{(s-s_5)^3}$，$s_1$有3个重根，$s_5$有2个重根。 \n", "令：$X_1(s)=\\frac{U(s)}{(s-s_1)^3}=\\frac{X_2(s)}{s-s_1}\\Rightarrow\\dot X_1=s_1X_1+S_2$\n", "$$\\Rightarrow \\begin{cases}\\dot X=\\begin{bmatrix} s_1&1\\\\&s_1&1\\\\&&s_1\\\\&&&s_4\\\\&&&&s_5&1\\\\&&&&&s_5\\end{bmatrix}X+\\begin{bmatrix}0\\\\0\\\\1\\\\1\\\\0\\\\1\\end{bmatrix}u \\\\ Y=\\begin{bmatrix}k_{11}&k_{12}&\\cdots&k_{52}\\end{bmatrix}X\\end{cases}$$\n", "- 如果存在**零极点对（零极点重合）**，则相应零极点对相消。 \n", "``` matlab\n", "sys_min=minreal(sys)\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**显然，我们可以发现，对于输入量含有导数项的物理方程，我们通过构建输入输出传递函数获得状态空间模型更加容易一些。**\n", "\n", "**但是，输入输出传递函数模型只能建立单输入单输出系统(SISO)，因此，通常对于多输入多输出系统(MIMO)，我们需要对每一个输入输出的组合构建一个传递函数模型**\n", "\n", "### Example\n", "\n", "仍对于这个高阶微分方程$\\dddot y + 3\\ddot y + 5\\dot y + 7y = 2\\dot u + 4u$，拉普拉斯变换后，得：\n", "$$G(s)=\\frac{Y(s)}{U(s)}=\\frac{2s+4}{s^3+3s^2+5s+7}$$\n", "有：\n", "``` matlab\n", "num=[2 4];den=[1 3 5 7];\n", "```" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "system_1 =\n", " \n", " 2 s + 4\n", " ---------------------\n", " s^3 + 3 s^2 + 5 s + 7\n", " \n", "Continuous-time transfer function.\n", "\n", "\n" ] } ], "source": [ "num=[2 4];den=[1 3 5 7];\n", "system_1=tf(num,den)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 状态空间模型转化为传递函数模型\n", "\n", "状态空间模型：$$\\begin{cases}\\dot X = AX + BU \\\\Y = CX + DU \\end{cases}\\Rightarrow \\begin{cases}sX = AX + BU \\\\Y = CX + DU \\end{cases}$$ \n", "传递函数模型：$$G(s)=\\frac{Y(s)}{U(s)}=C(sI-A)^{-1}B+D$$ " ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "system_2 =\n", " \n", " 2 s + 4\n", " ---------------------\n", " s^3 + 3 s^2 + 5 s + 7\n", " \n", "Continuous-time transfer function.\n", "\n", "\n" ] } ], "source": [ "system_2=tf(sys)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## - 流程图模型\n", "\n", "以双输入双输出为例： \n", "$\\begin{cases}\\ddot y_1 +a_1\\dot y_1 +a_2y_2 =b_1\\dot u_1 +b_2u_1 +b_3u_2 \\\\ \\dot y_2 +a_3y_2 +a_4y_1 =b_4u_2\\end{cases}$ \n", "$\\Rightarrow$\n", "$\\begin{cases}y_1=\\lmoustache(-a_1y_1+b_1u_1)dt +\\lmoustache\\lmoustache(b_2u_1+b_3u_2-a_2y_2)dt^2 \\\\ y_2=\\lmoustache(-a_3y_2-a_4y_1+b_4u_2)dt\\end{cases}$\n", "\n", "![流程图]() \n", "把每一次积分结果看作是状态变量得一项，则： \n", "$\\begin{cases}\\dot x_1=-a_1x_1+x_2+b_1u_1 \\\\ \\dot x_2=-a_2x_3+u_1b_2+u_2b_3 \\\\ \\dot x_3=-a_4x_1-a_3x_3+u_2b_4 \\\\ y_1=x_1\\\\y_2=x_3\\end{cases}$ \n", "$\\Rightarrow$\n", "$\\begin{cases}\\dot X=\\begin{bmatrix}-a_1&1&0\\\\0&0&-a_2\\\\-a_4&0&-a_3\\end{bmatrix}X +\\begin{bmatrix}b_1&0\\\\b_2&b_3\\\\0&b_4\\end{bmatrix} \\begin{pmatrix}u_1\\\\u_2\\end{pmatrix} \\\\ Y=\\begin{bmatrix}1&0&0\\\\0&0&1\\end{bmatrix}X\\end{cases}$ \n", "\n", "**显然，对于MIMO系统，相比于输入输出传递函数模型，流程图模型，更容易构建状态空间模型。**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 线性系统组合模型\n", "\n", "- 并联：$G(s)=G_1(s)+G_2(s)$ \n", "```matlab\n", "parallel(sys1,sys2) \n", "```\n", "- 串联：$G(s)=G_1(s)G_2(s)$ \n", "```matlab\n", "series(sys1,sys2)\n", "```\n", "- 反馈：$G(s)=(I+G_1(s)G_2(s))^{-1}G_1(s)$ \n", "```matlab\n", "feedback(sys1,sys2)\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## - 非线性系统模型\n", "\n", "相比于线性系统，非线性系统是更加普遍的系统。通常非线性系统是无法直接地写出其状态空间模型。它首先需要被线性化，我们才能在某一个平衡点附近写出近似状态空间模型。 \n", "以动力学方程为例：\n", "$$A(q)\\ddot q+C(q,\\dot q)\\dot q+g(q)=Qu$$\n", "状态变量$x=[q,\\dot q]^T$在平衡点$x_e=[q_e,0]^T$处有稳定输入 $u_e$。\n", "则：\n", "$$\\dot X=F(x)+G(x)u=\\begin{bmatrix} \\dot q \\\\ A(q)^{-1}(-C(q,\\dot q)\\dot q-g(q))\\end{bmatrix} +\\begin{bmatrix} 0 \\\\ A(q)^{-1}Q\\end{bmatrix}u$$\n", "\n", "$$\\Rightarrow\n", "\\begin{cases}A=\\frac{\\partial F}{\\partial q}(x_e)+\\frac{\\partial G}{\\partial q}(x_e)u_e=\\begin{bmatrix}O&I\\\\-A^{-1}\\frac{\\partial G}{\\partial q}({x_e})&O\\end{bmatrix}\\\\ B=G(x_e)=\\begin{bmatrix}O\\\\A^{-1}Q(x_e)\\end{bmatrix}\\\\g(x_e)=Qu_e\\end{cases}$$\n", "\n", "### Example 倒立摆\n", "$$\\begin{cases} (J+ml^2)\\frac{d^2\\theta}{dt^2}+(mlcos\\theta)\\frac{d^2x}{dt^2}=mldsin\\theta \\\\ (M+m)\\frac{d^2x}{dt^2}+\\mu_c\\frac{dx}{dt}+(mlcos\\theta)\\frac{d^2\\theta}{dt^2}=u+mlsin\\theta(\\frac{d\\theta}{dt})^2\\end{cases}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 状态空间模型的线性变换和Jordan规范型\n", "\n", "为了降低系统分析和设计的难度，把一般形式的状态空间模型变换成为特定形式的模型，并且不改变系统本身的特性（传输特性），这个过程被称作状态空间模型的线性变换。Jordan规范型是其中一种特定形式的模型。\n", "\n", "- **线性变换** \n", "原始状态空间模型：$$\\sum(A,B,C,D):\\begin{cases}\\dot X = AX + BU \\\\Y = CX + DU \\end{cases}$$\n", "经线性变换：$\\bar{X}=PX$ \n", "得：$$\\sum(\\bar{A},\\bar{B},\\bar{C},\\bar{D}):\\begin{cases} \\dot {\\bar{X}} = \\bar{A} \\bar{X} + \\bar{B} U \\\\ Y = \\bar{C} \\bar{X} + \\bar{D} U \\end{cases}$$\n", "则：\n", "$$\\bar{A}=PAP^{-1},\\bar{B}=PB,\\bar{C}=CP^{-1},bar{D}=D$$\n", "那么：$$\\bar{G}=\\frac{Y}{U}=\\bar{C}(sI-\\bar{A})^{-1}\\bar{B}+\\bar{D}=C(sI-A)^{-1}B+D=G$$\n", "matlab函数为：\n", "``` matlab\n", "sys_new=ss2ss(sys_old,P)\n", "```\n", "- **对角线规范化** \n", "将系统矩阵A化为对角线矩阵，意味着$\\bar{A}$为矩阵A得特征值矩阵，因此，利用matlab函数：\n", "``` matlab\n", "[P,D]=eig(A)\n", "```\n", "得到对角线矩阵$D=\\bar{A}=PAP^{-1}$ \n", "- **Jordan规范化** \n", "由matlab函数：\n", "``` matlab\n", "[P,J]=jordan(A)\n", "```\n", "得到Jordan矩阵$J=\\bar{A}=PAP^{-1}$" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "A =\n", "\n", " 1 -3 -2\n", " -1 1 -1\n", " 2 4 5\n", "\n", "\n", "P =\n", "\n", " -1 1 -1\n", " -1 0 0\n", " 2 0 1\n", "\n", "\n", "J =\n", "\n", " 2 1 0\n", " 0 2 0\n", " 0 0 3\n", "\n", "\n" ] } ], "source": [ "A = [1 -3 -2; -1 1 -1; 2 4 5]\n", "[P,J]=jordan(A)" ] } ], "metadata": { "kernelspec": { "display_name": "Matlab", "language": "matlab", "name": "matlab" } }, "nbformat": 4, "nbformat_minor": 2 }