微分方程（Differential Equation）是描述一个动态系统的时域数学模型，我们借此来间接研究系统的行为。微分方程的性质和解结构，与系统的响应特点有本质上的联系，这点需要在进行以传递函数为主的频域模型讲解前说明。

本文探究三个方面，ODE与LTI系统响应的关系，以及线性化的方法与合理性。后者是本文的重点内容。文中多次出现系统响应，系统输出，解轨迹都是指对应微分方程的解。

微分方程有很多种，我们在控制理论中多见的是常微分方程（Ordinary Differential Equation）。我们可以通过机理分析或者系统辨识来得到ODEs，从而进行基于模型的控制设计。

ODE分为线性(Linear)和非线性(Nonlinear)的。在经典控制的范畴中，我们关注点主要在线性常系数ODE，也会设计到非常系数的ODE。前者描述的系统我们称之为线性时不变系统，或者常称为线性定常系统（Linear Time-Invariant System，LTI System）。下面所述的ODE均为线性常系数ODE（Constant-Coefficient Linear ODE）

既然不是所有系统都是线性的，或者说，现实中的系统几乎都存在非线性。利用线性常系数ODEs来描述系统总是有一定误差的。因此我们自然想知道使用线性常系数ODEs的合理性。

ODE与系统响应

ODEs描述了系统的运动规律，ODE的解对应的就是某个系统的响应。我们所熟知的牛顿第二运动定律（Newton's second law of motion）就是由一个微分方程来刻画的

这是一个十分简单的二阶ODE，描述力物体受力后的位移变化规律。我们可以直接通过积分来得到其解。不过从一般解微分方程的角度来看，此方程的解依旧由两部分组成： 通解（general solution）以及特解（particular solution）。

特解： 目测可得其中一个特解 

由此可知全解为 

给定初始条件  ，那么

此微分方程的解为： 

特别的，当零初始条件时，解也就是高中物理里常用的公式 

那么当初始状态为0时，我们称这样的解是被该ODE描述的系统的零状态响应（Zero State Response ），此时系统的响应变化完全是由输入F所驱动的，所以也可以理解为纯输入响应。 零状态响应只需要将零初始条件赋予系统全解形式即可得到，本质上还是通解与特解之和，这里的例子比较特殊，通解部分正好为0。

现在在某一时刻，我们突然撤走F，系统该如何变化呢？ 此时描述系统的ODE的右端F变成了0，由非齐次ODE变成了齐次ODE

我们知道此时齐次系统的全解，也就是原非齐次的通解。假设初始条件为

那么对应的解为

我们称这样的解为该ODE描述的系统的零输入响应（Zero Input Response）。从撤走外力的那一刻起，系统的响应完全由系统自身的状态决定。试想如果此时系统的初始状态也为0，那么x就会从一开始就保持为0。

ODE解的叠加性

和  的零初始条件解  相加

得到的解：

对比(1)(2) 完全一致。从这个例子不足以推广到所有线性ODE，但是还是不妨碍我们搬出这个结论：

非齐次线性ODE的解，由其零初始条件下的非齐次方程的解 和 非零初始条件下齐次方程的解组成。其对应的系统的全响应，则由这两部分对应的零状态响应与零输入响应相加而得。

线性常系数ODE 全解 → 齐次方程非零初始条件解 + 非齐次方程零初始条件解

LTI系统Complete response → Zero State Response+ Zero Input Response

对于其他线性常系数ODE来说，这一叠加性也都是满足的。系统的解与ODE之间的关系，我们在下一章的传递函数和稳定性中会用到，对于理解传递函数的一些概念尤其重要。

ODE线性化方法

我们会遇到一些非线性的ODE，比如加入了一些三角函数或高阶多项式

虽然这仍然是一个二阶ODE，但由于是非线性的，并不满足叠加原理。对于这个方程描述的非线性系统，我们仍然可以通过线性化的方法来找到一个它的线性版本。

线性化的理论基础是泰勒展开（Taylor expansion），任意函数在定义域内某个点可以写成多项式的求和形式。于是我们只取泰勒展开的一阶导数项和常数项，作为原函数的近似。比如取展开点为x=0，上面式子中的cos(x)和x3根据泰勒公式有：

略去高次多项式，只取其线性部分，我们有如下近似

这样上述非线性系统就近似成了x=0附近的线性系统

ODE线性化的合理性

线性系统的这种线性化（Linearization）近似合不合理呢？ 从结果上来看是合理的，线性理论造就了一大批控制系统，连飞机都用到了线性化模型，实践看来是ok的。那么我们自然要问，这种近似是任何时候都成立的吗？ 答案是，No。

数学上我们知道只要x取值在展开点附近，泰勒展开的线性部分对原函数的近似还是令人满意的。从上节线性化的方法中可以看，展开点似乎是可以任意选择的，理论上可以把一个非线性ODE在不同的展开点处近似成无数个线性ODE。但线性方程与原非线性方程之间的近似关系只有在系统的解轨迹（trajectory）不远离展开点附近的情况下才成立。根据泰勒公式我们知道，当  ，即系统某时刻的解与展开点差的模，为一个比较小的数时，它的高阶多项式才能忽略，线性化的近似关系才能成立。否则，高次项的值远大于线性部分，线性方程描述的系统与原非线性方程描述的系统已经相去甚远。

这样产生了两个问题：

1. 原非线性方程解空间中这么多点，究竟要在哪个点进行线性化？
2. 线性化后得到的线性方程，怎么保证解就能在原来展开点的附近呢？

我们尝试回答这两个问题

一，哪个点进行线性化，取决于我们希望系统稳态时所处的工作点（set point）。

在regulation问题中我们希望系统的输出能够稳定到set point。一旦有了外来扰动，那么控制器就会给出控制信号来继续稳定输出。控制器设计指标之一就是不让输出有大幅度变化，要又快，又稳，又准地把输出再维持到设定值。既然regulation中我们输出一旦稳定便不再大幅度变化，我们就可以把线性化的展开点放到这个set point上，这个稳态时的点，也叫做系统的稳定平衡点（stable equilibrium）。具体可以这么做：

1. 选择一个期望的系统稳定时的set point
2. 对系统的非线性ODE在该工作点处进行线性化。
3. 按照该线性系统的模型设计控制器，控制非线性系统的输出（output），或者也说响应（response）稳定到这个工作点。

带入原方程得到

同样选择让误差e稳定到你希望的set point，一般都是0，来线性化这个非线性ODE。接下来的工作就是设计控制器。

二，控制器的设计以及扰动对线性化模型的有效性起了很大影响。

我们注意到上小节中的第3步，按照线性化模型设计控制器时，我们就假定了线性模型和原系统是可以近似等价的，但等价是有条件的。反过来设计出来的控制器，如果无法满足近似等价的条件，那这时候控制器就可能要失效了，或者说是因为线性化模型已经不能很好地描述原来的非线性系统了。

我们说，选择了某个set point进行线性化后，线性化模型与原系统之间拥有良好近似关系的关键就在于控制器能否把系统的输出都控制在set point附近。另一方面，控制器controller的一大任务就是抗干扰。如果干扰特别大，可以想象输出很容易大幅度偏离set point，也会威胁线性化模型效果。

控制器如何保证受扰后的系统依旧能够把解稳定在set point附近呢？ 我们先说一个不加控制作用的系统，它天然地也可以具有把解稳定到自身稳定平衡点的作用。

来看非线性系统的例子，不加任何外力F的系统，则变成了齐次方程：

“这蜜汁漩涡…” 好像很密很乱，但是我们可以看到有一些解的轨迹最后“转转转”，落入了原点（origin）。而明显，我只取了x1 为0到4，那外围那些大圈圈是个什么鬼？

于是我让x（0）从3变化到4，步长取0.5，只有三条轨迹，这下总清晰了吧。

红色部分就是吸引域的估计，只要初值落在红色区域内 [Note2] ，那么方程的解最后就会落入原点origin。在外面的初值会导致解发散到无穷远。

那么我们能否改变系统的稳定平衡点，来满足我们的需求呢？

这个例子让我们意识到，设计控制器的目的其实就是改变系统结构，控制器的任务之一就是把系统的稳定平衡点给“挪”到我们自己设定的set point上。如果系统本来就没有这样的稳定平衡点，控制器的任务之一就是“创造”这样的稳定平衡点  [Note3]  。加了控制器的系统能够拥有上述非线性系统的特性，即把set point作为系统的稳定平衡点，把set point附近的解通通吸引过去。这样就能做到解总是在set point附近，从而保证了按照线性系统设计的控制器对非线性系统能够有效地控制。

线性控制系统的响应一般都会控制在小范围的波动。一旦由于外界干扰，或者参考输入信号突然的大幅度改变，导致误差骤增，就容易导致系统的不稳定，这十分考验控制器的能力。这是因为一旦误差骤增，与控制器自身的增益相乘后，将会得到一个数值很大的控制信号。这样的信号是为了让系统能够迅速修正误差，重新进入稳态。这点很容易理解，就好像你猛踩油门，为了达到120码一样。但这可能造成系统输出在短时间内的大幅度的波动，从而让解在某一时刻跑出被控系统的吸引域。一旦跑出吸引域，麻烦就大了，这时候解可能就再也跑不回set point了，而解持续增大，从而线性化模型就失去了近似的意义，由线性化模型设计的控制器就很有可能会失效。

实际系统的输入突然变得很大，不管何种类型的输入：控制器增益K太大导致控制器输出的信号过大，还是外界突然来了一个强干扰绕使得干扰输入过大，又或者一个参考信号突然丢失从而误差变大，这都会导致plant或者process的输入骤然变大，从而使得输出跑出其稳定平衡点的吸引范围。即便这些问题修复了，系统从那一时刻起的初始条件已然不会再让系统的输出返回稳定平衡点，也就是set point了。系统的非线性特性可能随着这种输入的突然增大而被激发，很有可能会让线性控制器失效。对于更复杂的非线性系统，情况会更加复杂，它们可能会拥有多个平衡点，我们不希望解最后落入我们不想让它去的地方。这些日后有机会再提。

间接地我们得到一个信息，控制器的一大设计目标，自然就是扩大吸引域的范围了。换句话说，这也就是控制系统鲁棒性（robustness）的表现，因为外界激励的增大导致线性模型失效其实也等同于模型参数的改变。如果不能实时地更新线性模型的参数和控制器的参数，那么这样的控制器就不能适应较大模型变动。好点的情况是控制器依旧能够稳定系统，但是控制性能指标得不到满足，差点的情况就是直接无法稳定而导致系统fail。

不光光是线性控制器，即便是非线性控制器，也需要尽量扩大吸引域。很多时候很难通过设计来赋予系统一个全局稳定平衡点（globally stable equilibrium） [Note4] ，即无论从哪个初值出发都最后落入稳定平衡点。这些就涉及到非线性控制（nonlinear control）和非线性系统（nonlinear system）上去了，这里不再多展开了。

总结一下

线性常系数ODE 的解描述了LTI系统的响应，他们之间有对应关系。

系统的线性化理论基础来自泰勒展开和小误差近似。

非线性系统的近似线性化，其初值条件一般不能偏离展开点太远，否则容易产生大偏差，而导致线性化模型不再准确，导致线性控制器的控制效果变差。甚至会激发原系统的非线性特性，控制器失效，导致系统变得不稳定。较大的干扰往往容易引发线性控制器的失效，此时引入非线性控制器可以较好地控制干扰，帮助系统的工作区域能扩展到更大的范围而保持稳定。

Note：

[1] 这个名字取得并不草率，而是相当专业的非线性动力学名词。

[2] 严格地说不能随便下这种结论，但是对于这里例子平衡点确实只有一个，所以这些红点大致覆盖了吸引域。对于有多个平衡点的非线性系统，粗略地画图是不能判断这片红色区域里还有没有其他不稳定平衡点的。

[3] 这里涉及到了一个是不是任何系统都能这样挪的问题，要考虑能控性，这里先不管。

[4] 线性系统本身只有一个平衡点，确实是全局稳定的平衡点，即吸引域无穷大。我们之后会再提。

Reference

[1] Hassan K.Khali，Nonlinear Systems , 3rd Edition, 2014, Pearson

[2] JJ E.Slotine, Weiping Li, Applied Nonlinear Control, 1991, Prentice Hall

[3] 胡寿松，自动控制原理（第六版），2013，科学出版社

Appendix

MATLAB codes for estimation of region of attraction(ROA)

% SJM Frank
% estimation of ROA
% 06-22-2019
dt = 0.01;        %step size
t_end = 50;
t = 0:dt:t_end;   %time sequence
N = size(t,2);

y =zeros(2,N);
for j = -10:0.2:10         % initial value for x1
    for k = -10:0.2:10     % initial value for x2
        x =zeros(2,N);
        x(:,1)=[j;k];
   
        for i = 1 : N-1     %Euler method
            x(1,i+1) = x(1,i) + x(2,i) * dt;
            x(2,i+1) = x(2,i)+ (-cos(x(1,i)).*x(2,i)-x(1,i)-x(1,i).^3) * dt;
        end
    hold on
    if abs (x(1,end))<0.1 ||abs( x(2,end))<0.1  % if converges. 0.1 is changable.
        plot(j,k,'xr')
    end
    end
end
ylim([-10 10])
xlim([-10 10])
xlabel('x1')
ylabel('x2')
% legend('x(0)=3','x(0)=3.5','x(0)=4')
% figure;
% plot(y(1,:),y(2,:))

若有纰漏，烦请指出。转载还请私信联系，请勿私自转载。O(∩_∩)O