状态空间

ABCD形式

回顾前面关于常微分方程的讨论,我们可将微分方程表述为如下形式:

\[\begin{split}\dot{\vec{x}}(t) &= \vec{f}(\vec{x}(t), \vec{u}(t), t) \\ \vec{y}(t) &= \vec{g}(\vec{x}(t), \vec{u}(t), t)\end{split}\]

在这种形式下,\(x\)表示系统中的状态,\(u\)表示系统的任何外部指定的输入,而\(y\)表示该系统的输出(亦即可能并非状态的变量,但这些变量可通过状态和输入的值求得。)

这些方程有一个有趣的特例。这个特例发生在是当函数\(\vec{f}\)以及\(\vec{g}\)和向量\(\vec{x}\)还有\(\vec{u}\)线性相关时。在这种情况下,方程可以改写为:

\[\begin{split}\dot{\vec{x}}(t) &= A(t) \vec{x}(t) + B(t) \vec{u}(t) \\ \vec{y}(t) &= C(t) \vec{x}(t) + D(t) \vec{u}(t)\end{split}\]

此问题中的矩阵是所谓“ABCD”矩阵。“ABCD”形式是很有用的,因为如果系统用这种形式表达的,那么我们就可以在系统上进行一些有趣的计算。例如,我们可以使用\(A\)矩阵计算系统的自然频率。而使用这些矩阵的不同组合,我们可以确定对系统控制相关的几个极其重要的性质(例如可观察性和可控性)。

请注意,ABCD形式允许这些矩阵随时间变化。有一种稍为更专门的形式,除去线性以外,也不随时间变化:

\[\begin{split}\dot{\vec{x}}(t) &= A \vec{x}(t) + B \vec{u}(t) \\ \vec{y}(t) &= C \vec{x}(t) + D \vec{u}(t)\end{split}\]

这种形式通常被称LTI(线性时不变)形式。LTI形式是很重要的。因为这种形式除了具有与ABCD的形式一样的特殊性质外, 还可以以一个非常简单的形式实现“模型交换”。以前,当用户(采用手工方法或者用建模软件)推导给定系统的行为方程后,其中一个将这些方程导入其他工具的方法,就是把方程变为LTI形式。这意味着,该模型可以用带有数字或表达式的一系列矩阵进行交换、共享或者发布。如今,新技术如Modelica语言和FMI为模型交换提供了更好的选择。

LTI模型

如果有人提供了我们一个LTI形式的模型,我们要如何用Modelica描述它呢?这是其中一种可能的方法:

model LTI
  "Equations written in ABCD form where matrices are also time-invariant"
  parameter Integer nx=0 "Number of states";
  parameter Integer nu=0 "Number of inputs";
  parameter Integer ny=0 "Number of outputs";
  parameter Real A[nx,nx]=fill(0,nx,nx);
  parameter Real B[nx,nu]=fill(0,nx,nu);
  parameter Real C[ny,nx]=fill(0,ny,nx);
  parameter Real D[ny,nu]=fill(0,ny,nu);
  parameter Real x0[nx]=fill(0,nx) "Initial conditions";
  Real x[nx] "State vector";
  Real u[nu] "Input vector";
  Real y[ny] "Output vector";
initial equation
  x = x0 "Specify initial conditions";
equation
  der(x) = A*x+B*u;
  y = C*x+D*u;
end LTI;

第一步,是在模型里声明nxnuny等参数。这些参数分别代表了状态、输入和输出的数量。然后,我们定义矩阵ABC以及D。因为我们正在创建一个线性时不变表示下的模型,所以所有这些矩阵都可以被定义为参数。由于矩阵ABCD的定义紧接着 [ 以及 ] ,因此我们知道它们均是数组。而由于[]内有标示了两个维度,我们知道以上的数组均是矩阵。最后,我们看到x0xu以及y的变量声明。这些变量也都是数组。不过,由于这些变量仅有一维,因此它们都是向量。

本模型的另一个特点是其参数均有默认值。对于nxnuny等参数,默认的假设是状态、输入和输出的数量均为零。而对于矩阵则默认其元素均为零。除非另有规定,对于初始条件我们同样假设所有的状态在仿真开始取零值。我们将很快看到,这些假设如何可以让我们通过直接覆盖参数值来写编写简单的模型。

向量方程

该模型的其余部分现在看起来应该非常熟悉了。必须指出,模型中的方程皆是矢量方程。Modelica语言中方程可以包括标量或数组。对于矢量方程唯一的要求是,方程的两侧需要有相同的维数和以及每个维度大小均相等。因此,在本例的LTI模型里,我们有以下的初始化方程:

initial equation
  x = x0 "Specify initial conditions";

这个方程是一个向量方程,内容是x的每一个元素在仿真开始时等于其在x0的对应元素值。实际上,这些向量中的每组对应元素会自动展开为一系列的标量方程。

还有另外一点有助于保持这些方程可读性。Modelica语言包含了关于函数向量化的一些特殊的规则。概括地说,这些规则规定了,倘若你有一个函数可以对标量进行运算,那么你就可以立刻用这个函数可以进行向量运算。如果你尝试用该函数进行向量运算, Modelica语言会自动将函数应用在向量中的每个元素上。因此,LTI模型内的der(x)表达式就是一个表示x中每个元素微分的向量。

最后,许多代数运算符(如:+-以及*)在应用于向量或矩阵时有着特殊的意义。运算符的定义设计得与常规的数学符号对应。所以,在LTI模型里,表达式A*x对应了矩阵与向量的积。

LTI例子

考虑到这一点的所有,让我们重新审视我们几个以前的例子,看看他们如何能在长期激励形式表示使用我们的LTI模型。讨论了上述的内容后,让我们重新检视前述的数个例子。目的是看看这些可以如何使用LTI模型表示为LTI形式。注意,我们会再次使用继承(通过使用extends关键词)以重用LTI模型的代码。

让我们从前面介绍的简单的一阶系统开始。使用LTI模型,我们可以将模型重写为:

model FirstOrder "Represent der(x) = 1-x"
  extends LTI(nx=1,nu=1,A=[-1], B=[1]);
equation
  u = {1};
end FirstOrder;

在我们扩展继承LTI时,我们只需要指定与默认值不同的参数即可。在这里,我们指定模型有一个状态和一个输入。然后,我们定义AB为1x1矩阵。最后,由于有一个输入,我们需要为这个输入提供一个方程。该输入一般而言可以随时间变化,因此我们不把它表示为参数,而表示为方程。请注意,在方程中有:

u = {1};

表达式{1}是一个向量源代码文本。这表示,我们其元素组成的列表来构建向量。在这里,向量仅有一个元素1。但我们建立一个用逗号分隔开的表达式列表去创建较长的向量,例如:

v = {1, 2, 3*4, 5*sin(time)};

值得一提的是,我们在extends语句内除了设置参数值,也可以包含等式。因此,我们可以完全避免equation区域,而将模型简化为:

model FirstOrder_Compact "Represent der(x) = 1-x"
  extends LTI(nx=1,nu=1,A=[-1], B=[1], u={1});
end FirstOrder_Compact;

一般来说,加入equation区域可以使代码有点更具可读性。但也有一些情况下,向extends语句加入等式作为对模型的修改会更为方便。

现在,让我们来关注也是前面讨论过的冷却模型。我们可以把模型用LTI形式改写如下:

model NewtonCooling "NewtonCooling model in state space form"
  parameter Real T_inf=27.5 "Ambient temperature";
  parameter Real T0=20 "Initial temperature";
  parameter Real hA=0.7 "Convective cooling coefficient * area";
  parameter Real m=0.1 "Mass of thermal capacitance";
  parameter Real c_p=1.2 "Specific heat";
  extends LTI(nx=1,nu=1,A=[-hA/(m*c_p)],B=[hA/(m*c_p)],x0={20});
equation
  u = {T_inf};
end NewtonCooling;

这个模型非常类似于前一个。然而,在这种情况下我们并不是把数字直接输入矩阵里。相反,我们用输入带有参数mc_p等等的表达式。这样的话,当这些物理参数改变时,矩阵AB的值也会相应改变。

我们可以采取类似的做法把前面的机械示例改写为LTI形式。

model RotationalSMD
  "State space version of a rotational spring-mass-damper system"
  parameter Real J1=0.4;
  parameter Real J2=1.0;
  parameter Real k1=11;
  parameter Real k2=5;
  parameter Real d1=0.2;
  parameter Real d2=1.0;
  extends LTI(nx=4, nu=0, ny=0, x0={0, 1, 0, 0},
                  A=[0, 0, 1, 0;
                     0, 0, 0, 1;
                     -k1/J1, k1/J1, -d1/J1, d1/J1;
                     k1/J2, -k1/J2-k2/J2, d1/J2, -d1/J2-d2/J2]);
equation
  u = fill(0, 0);
end RotationalSMD;

同样,我们从物理参数得到A的值。在本例里要注意A的构造。数学上,矩阵\(A\)被定义为:

\[\begin{split}A &= \left| \begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -\frac{k_1}{J_1} & \frac{k_1}{J_1} & -\frac{d_1}{J_1} & \frac{d_1}{J_1} \\ \frac{k_1}{J_2} & -\frac{k_1}{J_2}-\frac{k_2}{J_2} & \frac{d_1}{J_2} & -\frac{d_1}{J_2}-\frac{d_2}{J_2} \\ \end{array} \right|\end{split}\]

在构造\(A\)时,我们可以注意到其前两列可以更容易表示为一个零矩阵以及一个单位矩阵的组合。换句话说,将矩阵表示为子矩阵的组合可能更为清晰,即:

\[\begin{split}A &= \left| \begin{array}{cc} \left| \begin{array}{cc} 0 & 0 \\ 0 & 0 \end{array} \right| ~ \left| \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right| \\ \left| \begin{array}{cc} -\frac{k_1}{J_1} & \frac{k_1}{J_1} \\ \frac{k_1}{J_2} & -\frac{k_1}{J_2}-\frac{k_2}{J_2} \end{array} \right| \left| \begin{array}{cc} -\frac{d_1}{J_1} & \frac{d_1}{J_1} \\ \frac{d_1}{J_2} & -\frac{d_1}{J_2}-\frac{d_2}{J_2} \end{array} \right| \end{array} \right|\end{split}\]

在Modelica语言,我们可以如下地用子矩阵构建A矩阵:

model RotationalSMD_Concat
  "State space version of a rotationals spring-mass-damper system using concatenation"
  parameter Real J1=0.4;
  parameter Real J2=1.0;
  parameter Real k1=11;
  parameter Real k2=5;
  parameter Real d1=0.2;
  parameter Real d2=1.0;
  parameter Real S[2,2] = [-1/J1, 1/J1; 1/J2, -1/J2];
  extends LTI(nx=4, nu=0, ny=0, x0={0, 1, 0, 0},
                  A=[zeros(2, 2), identity(2);
                     k1*S+[0,0;0,-k2/J2], d1*S+[0,0;0,-d2/J2]],
                  B=fill(0, 4, 0), C=fill(0, 0, 4),
                  D=fill(0, 0, 0));
equation
  u = fill(0, 0);
end RotationalSMD_Concat;

本节里,我们并没有把Lotka-Volterra方程表示为LTI形式。Lotka-Volterra方程虽然是时不变系统,但是同时却是非线性的。值得一提的是,Modelica并没有在LTI模型内执行对线性和时不变这两个属性的约束。因此,用上述方法其实可以描述非线性或者时变的模型。但倘若如此就会引起一定的混乱,因为LTI这个术语意味这方程是线性时不变的。

使用部件

在至今所有的例子中,我们已经(通过extends)用继承来重用LTI模型的公式。一般而言,把方程作为子组件是一个相比之下好得多的代码复用方法。为了说明这种方法,我们将前面讨论过的electrical examples重写为LTI形式。不过这次,我们会建立为LTI模型创建一个命名实例。

model RLC "State space version of an RLC circuit"
  parameter Real Vb=24;
  parameter Real L=1;
  parameter Real R=100;
  parameter Real C=1e-3;
  LTI rlc_comp(nx=2, nu=1, ny=2, x0={0,0},
               A=[-1/(R*C), 1/C; -1/L, 0],
               B=[0; 1/L],
               C=[1/R, 0; -1/R, 1],
               D=[0; 0]);
equation
  rlc_comp.u = {Vb};
end RLC;

请注意,这一次我们没有使用extends或任何形式的继承。相反,我们实际上声明了一个类型为LTI而名为rlc_comp的变量。一旦我们介绍完Modelica语言中描述不同行为的所有基础知识,我们就会将注意力转向如何将这些方程整理成可重用的组件。但现在,这不过是对后面的(重要)内容“先睹为快”而已。

我们在这个RLC例子中看到的是,我们现在有一个叫做rlc_comp的变量,而此部件拥有所有LTI模型的参数和变量。所以,例如我们可以看到用于指定输入u的方程写作:

rlc_comp.u = {Vb};

请注意,我们所提供的这个方程中变量urlc_comp里面的变量。正如我们将在后面看到的,我们可以用层次结构来管理在描述复杂系统时产生的复杂度。在这里,使用.操作符可以让我们引用层级结构里的变量。同样,我们会在介绍组件时对此进行彻底讨论。