第3章〓控制系统设计与仿真
建立控制系统数学模型的目的之一是对控制系统进行设计与分析，以期实现对工程对象良好的控制性能。本章结合运动控制与电力电子控制的8个具体工程问题，给出基于系统建模与仿真实验的系统设计与分析； 其中所涉及的理论与技术内容有助于我们深入理解已学过的知识，扩大自己的知识面。

3.1直流电动机转速/电流双闭环控制系统设计

自20世纪70年代以来，国内外在电气传动领域里，大量地采用了晶闸管整流电动机调速技术（简称VM调速系统）。尽管当今功率半导体变流技术已有了突飞猛进的发展，但在工业生产中VM系统的应用还是占有相当比重； 一般情况下，VM系统均设计成图31所示的转速/电流双闭环控制形式。



图31直流电动机双闭环调速系统结构图


1. 系统建模

根据2.5.5节的分析，可得图31所示的双闭环直流调速系统的动态系统结构，如图32所示。



图32双闭环调速系统的动态结构框图


2. 电流环与转速环调节器设计

（1） 双闭环控制的目的。

双闭环控制的直流调速系统着重解决了以下两方面的问题： 

① 启动的快速性问题。借助于PI调节器的饱和非线性特性，使得系统在电动机允许的过载能力下尽可能地快速启动。理想的电动机启动特性如图33所示。

②  提高系统抗扰性能。通过调节器的适当设计可使系统“转速环”对于电网电压及负载转矩的波动或突变等扰动予以控制（迅速抑制），在最大速降、恢复时间等指标上达到最佳，其动态特性如图34所示。




图33理想的电动机启动特性




图34双闭环控制直流调速系统负载扰动特性



（2） 关于积分调节器的饱和非线性问题。

双闭环VM调速系统中的ASR与ACR一般均采用PI调节器。在图35中给出了控制系统的PI控制规律动态过程，从中我们可知： 



图35比例积分调节器结构及其输入输出动态过程



① 只要偏差ΔU存在，调节器的输出控制电压Uc就会不断地无限制地增加。因此，必须在PI调节器输出端加限幅装置。

②  当ΔU=0时，Uc为常数。若要使Uc下降，必须使ΔU<0。因此，在直流调速控制系统中，若要使ASR退出饱和状态（进入线性控制状态），就一定要产生超调现象。

③  对于前向通道带有惯性环节的控制系统，若控制器存在积分作用，则在给定作用下，系统输出一定会出现超调。

（3） ASR与ACR的工程设计方法［1］。

对于直流电动机转速控制系统设计问题，通常应用“典型系统工程设计方法”； 对于电流环控制器的设计，在稳态上希望电流控制无静差，以得到理想的启动特性，同时要求电流的跟随性能要好； 因此，把电流环设计校正成典型Ⅰ型系统； 而转速环控制器，通常按把系统综合成典型Ⅱ型系统来设计，这样既可以保证转速无静差，又有较强的抗扰性。

对于图31所示的双闭环控制系统，理论上ASR与ACR均采用PI调节器的形式，且有如下最佳设计方法： 

① 电流调节器： 

WACR(s)=Kiτis+1τis

式中，Ki为电流调节器的比例系数； 
τi为电流调节器的超前时间常数。

为了让调节器零点对消掉控制对象的大时间常数（极点），选择

τi=Tl

在一般情况下，希望超调量σ%≤5%时，可取阻尼比ξ=0.707，KlTi=0.5，得

Kl=12Ti,Ti=Ts+Toi


又因为

Kl=KiKsβτiR


得到

Ki=KlτiRKsβ=TlR2KsβTi=0.5RKsβTlTi


②  转速调节器： 

WASR(s)=Knτns+1τns 

式中,Kn为转速调节器的比例系数； 
τn为转速调节器的超前时间常数。

转速开环增益

KN=KnαRτnβCeTm


按照典型Ⅱ型系统的参数选择方法


τn=hTn,Tn=2Ti+Ton
KN=h+12h2T2n

得到ASR的比例系数

Kn=(h+1)βCeTm2hαRTn


在工程上，通常选择h=5为最佳。
所以有

τn=5×Tn,KN=650×T2n

经过以上设计的VM调速系统，理论上讲有如下动态性能： 电动机启动过程中电流的超调量为4.3%，转速的超调量为8.3%； 需要说明的是，上述设计是以“线性系统”为前提条件的； 在系统调试中，受系统存在非线性、近似误差等因素的影响，实际动态性能会有所差异，下面的仿真实验也证明了这一点。

3. 仿真实验

（1） 仿真模型搭建。

系统中采用三相桥式晶闸管整流装置，基本参数如下： 

直流电动机： 220V，13.6A，1480r/min，Ce=0.131V/（r/min），允许过载倍数λ=1.5； 晶闸管整流装置： KS=76； 电枢回路总电阻： R=6.58Ω； 直流电动机时间常数： Tl=0.018s，Tm=0.25s； 反馈系数： α=0.00337V/（r/min），β=0.4V/A； 反馈滤波时间常数： Toi=0.005s，Ton=0.005s。

为使仿真结果更具真实性，系统控制对象的模型精度尤为重要； 这里，我们借助MATLAB/SimPowerSystems工具软件来建立系统的仿真模型。

根据图31的双闭环调速系统结构，应用SimPowerSystems中的晶闸管整流和触发装置、直流电动机模型作为系统的被控对象，可得晶闸管直流电动机调速系统如图36所示。



图36晶闸管直流电动机调速系统仿真结构图



① 转速环控制器/ASR。

由文献［1］知，转速环控制器需将外环系统校正成典型Ⅱ型系统，故转速环控制器应采用PI控制方式； 这里采用Simulink中的连续PID控制器模块，并且将微分项的系数设置为0，其结构如图37所示。



图37转速环调节器



图38所示为PID调节器内部参数的设置方法。



图38PID调节器参数设置


在PID Advanced选项卡中，限幅电压设为±8V（ASR输出限幅值=电机电流最大值×电流环反馈系数β值），抗积分饱和方法选择clamping法（在积分达到限幅值时停止积分），以降低系统的超调量，缩短调节时间。

注意，转速环控制器/ASR在系统的启动过程控制中，呈现“饱和非线性特性”，用以实现“最大电流/转矩启动”的最佳工作特性； 这一“非线性特性”在上面的理论设计中并未考虑。

② 电流环控制器/ACR。

由文献［1］知，电流环控制器需将内环系统校正成典型Ⅰ型系统，故电流环控制器也需采用PI控制方式，这里我们采用子系统的方式将电流环控制器封装为子系统。

需要强调的是，在SimPowerSystems环境下的电流环控制器与在Simulink环境下建立的数学模型不同，SimPowerSystems环境里的模型是物理模型，必须考虑实际情况，由于电流环的输出连接至晶闸管整流器的相位输入，所以需要将电流环的“电压输出转换为晶闸管的相位输入”进行控制。各电压值与各相位值为线性关系，故设线性方程为y=ax+b。

根据分析，负载电动机可以等效为阻感负载，当采用三相桥式全控整流电路时，相位控制范围为（0°~90°），对应的电压范围为（248.19V~0）。

另外，与采用线性模型仿真所设计的电流环调节器不同，将电压转换为对应的触发角度值，需要添加“限幅环节”，将角度值限幅在0°~90°之间。电流环子系统仿真程序如图39所示。



图39电流环子系统内部结构图


注意，通常电流环控制器/ACR也设置“饱和限幅输出特性”，以限定电机两端的最高电压； 但是，ACR应该始终工作在线性段，以保证电流控制的快速跟随特性。

③ 晶闸管触发和整流装置。

这里，我们采用三相桥式全控整流电路得到电动机驱动电压，三相桥式全控整流电路工作时要保证任何时候都有两只晶闸管导通，这样才能形成向负载供电的回路，并且是共阴极和共阳极组成各一个，不能为同一组的晶闸管； 因此，采用SimPowerSystems环境下的同步六脉冲发生器(Synchronized 6Pulse Generator)提供三相整流电路的触发脉冲。



图310晶闸管触发和整流装置图


图310中设置三个交流电压源Va、Vb、Vc相位角依次相差120°，得到整流桥的三相电源。采用电压测量模块测得线电压Vab、Vbc、Vca作为同步六脉冲发生器的输入端，alpha_deg为触发角输入端，用来接收电流环调节器的输出相位控制信号，同步六脉冲发生器输出相位是作用在6个晶闸管上的六脉冲向


图311电动机模块


量形式。为方便调试，我们采用多路测量电压表对各晶闸管两端电压及整流后的输出电压通过示波器进行观察。

④ 直流电动机。


本实验采用他励式励磁电动机，电动机模块如图311所示，该电动机的端子F+与F-分别接励磁电压源的正负极，A+与A-接输入电压，TL为电机负载输入，m端为输出端，电动机参数设置如图312所示。



图312电动机参数设置



⑤ 反馈滤波环节的滞后对消。

在实际工程中，为保证电流/电压反馈信号的质量，均采用“一阶惯性滤波器”（如图36中current feedback与speed feedback环节），为对消其产生的时间滞后，在前向通道上增加了Transfer Fcn1与Transfer Fcn2，以保证信号传输时间上的一致性。

（2） 系统仿真实验。

① 开环系统性能分析。

实验中，选择ode23tb或ode23t解算方法，并将系统powergui模块中Simulation type设置为Continuous，选取Start time=0.0，Stop time=6.0，仿真时间从0s到6s，晶闸管触发角初始为0°，1.6s时变为45°，3s时变为60°。

仿真模型如图313所示。



图313晶闸管直流电动机开环系统仿真模型



图314为晶闸管触发角给定和电动机转速动态特性仿真结果，通过仿真分析，我们可以看到电动机转速随晶闸管触发角的改变而变化的情况。



图314系统开环特性


上述仿真实验结果/实验曲线程序如下： 

程序1

clf

%数据读取

load speed.mat

t=signals(1,:);

y1=signals(2,:);

y2=signals(3,:);

%绘制曲线

［AX,H1,H2］=plotyy(t,y1,t,y2);

%坐标范围设置

set(AX(1),'xlim',［0,6］,'xTick',0:1:6);

set(AX(2),'xlim',［0,6］,'xTick',0:1:6);

set(AX(1),'ylim',［-10,90］,'yTick',-10:10:90);

set(AX(2),'ylim',［0,2000］,'yTick',0:400:2000);

%坐标轴颜色设置

set(AX(1),'XColor','k','YColor','k');

set(AX(2),'XColor','k','YColor','k');

%曲线颜色设置

set(H1,'color','b');

set(H2,'color','k');

%绘制网格

grid on;

%标注设置

set(get(AX(1),'Xlabel'),'string','时间(s)');

set(get(AX(1),'Ylabel'),'string','晶闸管出发角(°)');

set(get(AX(2),'Ylabel'),'string','电动机转速(r/min)');

title('晶闸管触发角与电动机转速');

gtext('晶闸管触发角');

gtext('电动机转速');

② 双闭环系统起动特性分析。

实验中，选择ode23tb或ode23t解算方法，并将系统Powergui模块中Simulation type设置为Discrete，采样周期设置为5e006s，最大仿真步长设置为5e006s，选取Start time=0.0，Stop time=3.0，仿真时间从0s到3s。

图315、图316、图317分别为ASR的输出与电动机转速动态特性仿真结果(其中图(b)为图(a)在0s附近的放大图)、ACR的输出与电动机转速动态特性仿真结果以及电动机电流与电动机转速动态特性仿真结果。



上述仿真实验结果/实验曲线程序参考前述的程序1。

通过仿真结果分析，我们可以看到，对于系统起动性能指标来说，起动过程中电流的超调量为5.7%，转速的超调量为0.2%； 这一结果与理论设计结果相近（或优于理论结果： 电流超调量为4.3%，转速超调量为8.3%），达到预期设计目的。

③ 双闭环系统抗扰性能分析。

实验中，我们选取Start time=0.0，Stop time=5.0，仿真时间从0s到5.0s。扰动加入的时间均为3.0s。

一般情况下，双闭环调速系统的干扰主要是负载突变与电网电压波动两种。图318绘出了该系统电动机转速在突加负载（ΔT=8N·m）情况下电动机电流Id与输出转速n的关系； 图319、图320分别绘出了电网电压突减（ΔU=100V）和电网电压突增（ΔU=100V）情况下晶闸管触发整流装置输出电压Ud0、电动机两端电压Ud，与输出转速n的关系。



上述仿真实验结果/实验曲线程序参考前述的程序1。

仿真结果分析如下。
对于该系统的抗扰性能，有如下几个结论： 

（1） 系统对负载的大幅度突变具有良好的抗扰能力，在ΔT=8N·m的情况下系统速降为Δn=12r/min（0.81%），恢复时间tf=0.25s（开环Δn=492r/min，tf=8.4s）。

（2） 系统对电网电压的大幅波动也同样具有良好的抗扰能力，在ΔU=100V的情况下，系统速降仅为6r/min（0.4%），恢复时间tf=0.3s。

（3） 该系统的起动时间和恢复时间（轻载情况下约为1.5s）都很短，基本能够满足工程上的设计要求。

可见双闭环控制系统较开环系统的性能有较大的提升，达到了预期设计指标。

注： ΔU是电网电压换算成整流后输出直流电压的等效值。

4. 小结

本节针对直流电动机的转速控制问题，给出了基于“典型系统工程设计方法”的双闭环PID控制方案，基于MATLAB/SimPowerSystems工具软件的仿真实验表明系统设计的有效性与可行性。同时，我们还注意到： 


图315ASR输出特性





图316ACR输出特性





图317电动机电流特性





图318突加负载抗扰特性





图319电网电压突加抗扰性能





图320电网电压突减抗扰性能


（1） 仿真实验结果与典型系统理论相近。
由图315可见，电流动态响应超调量为5.7%（理论值为4.3%），转速动态响应超调量为0.2%（理论值为8.3%）； 这是由于“典型系统工程设计方法”的理论分析与设计中作了一些简化/近似处理，未考虑转速控制器/ASR的“输出限幅特性”、整流装置的内阻Rrec值随温度的变化等问题，而是将系统近似为“线性定常系统”来设计。

（2） 受MATLAB/SimPowerSystems仿真工具软件模型库精度的影响，工程系统装置上的实际调试结果与仿真结果还将有一定误差。

总之，系统建模为控制系统的设计提供了依据，仿真工具软件为工程实现提供了保障，可有效地加快新技术与新产品的研发进程。

3.2一阶直线倒立摆系统的双闭环PID控制方案

在图321所示的一阶倒立摆控制系统中，通过检测小车位置与摆杆的摆动角，来适当控制驱动电机拖动力的大小，控制器由一台IPC完成。




图321一阶倒立摆控制系统



本节将借助于Simulink封装技术——子系统，在模型验证的基础上，采用双闭环PID控制方案，实现倒立摆位置伺服控制的数字仿真实验。

1. 系统建模

(1) 对象模型。

由2.5.1节式(234)知： 

① 一阶倒立摆精确模型为

x¨=(J+ml2)F+lm(J+ml2)sinθ·θ·2-m2l2gsinθcosθ(J+ml2)(m0+m)-m2l2cos2θ
θ¨=mlcosθ·F+m2l2sinθcosθ·θ·2-(m0+m)mlgsinθm2l2cos2θ-(m0+m)(J+ml2)


当小车的质量m0=1kg，倒摆振子的质量m=1kg，倒摆长度2l=0.6m，重力加速度取g=10m/s2时得

X¨=0.12F+0.036sinθ·θ·2-0.9sinθ·cosθ0.24-0.09cos2θθ¨=0.3cosθ·F+0.09sinθ·cosθ·θ·2-6sinθ0.09cos2θ-0.24

② 若只考虑θ在其工作点θ=0附近(-10°<θ<10°)的小范围变化，则可近似认为

θ·2≈0sinθ≈θcosθ≈1


由此得到简化的近似模型为

X¨=-6θ+0.8Fθ¨=40θ-2.0F

其等效动态结构图如图322所示。




图322一阶倒立摆系统动态结构图



(2) 电动机、驱动器及机械传动装置的模型。

假设选用日本松下电工MSMA021型小惯量交流伺服电动机，其有关参数如下： 


驱动电压： U=0～100V额定功率： PN=200W

额定转速： n=3000r/min       转动惯量： J=3×10-6kg·m2

额定转矩： TN=0.64N·m        最大转矩： TM=1.91N·m

电磁时间常数： Tl=0.001s       机电时间常数： Tm=0.003s


经传动机构变速后输出的拖动力F=0～16N； 与其配套的驱动器为MSDA021A1A，控制电压UDA=0～±10V。

忽略电机的空载转矩和系统摩擦，认为驱动器和机械传动装置均为纯比例环节，并假设这两个环节的增益分别为
Kd和Km。

对于交流伺服电动机，其传递函数可近似为

KvTmTls2+Tms+1

由于是小惯性电机,其时间常数Tl、Tm相对都很小，这样可以进一步将电动机模型近似等效为一个比例环节： Kv。

综上，电动机、驱动器、机械传动装置三个环节就可以成为一个比例环节： 

G(s)=KdKvKm=Ks

Ks=Fmax/Umax=16/10=1.6


2. 模型验证  

尽管上述数学模型系经机理建模得出，但其准确性(或正确性)还需运用一定的理论与方法加以验证，以保证以其为基础的仿真实验的有效性。模型验证的理论与方法是一专门技术，本书受篇幅所限不能深入阐述。下面给出的是一种必要条件法，即我们所进行的模型验证实验的结果是依据经验可以判定的，其正确的结果是正确的模型所应具备的必要性质。

(1) Simulink子系统。

如同大多数程序设计语言中的子程序(如： C语言中的函数、MATLAB中的M文件)功能，Simulink中也有类似的功能——子系统。子系统通过将大的复杂的模型分割成几个小的模型系统，使得整个系统模型更加简捷，可读性更强。

把已存在的 Simulink模型中的某个部分或全部“封装”成子系统的操作程序如下： 

① 首先使用范围框将要“封装”成子系统的部分选中，包括模块和信号线(如图323所示)。为了使范围框圈住所需要的模块，常常需要事先重新安排各模块的位置(注意： 这里只能用范围框，而不能用Shift逐个选定)。



图323选择要封装的模块


② 在模块窗口菜单选项中选择Edit→Create Subsystem，Simulink将会用一个子系统模块代替选中的模块组，如图324所示。



图324封装后的模型图


③ 所得子系统模块将有默认的输入端口和输出端口。输入端口和输出端口的默认名称分别为In1和Out1。调整子系统和模型窗口的大小使之更加美观，如图325所示。



图325调整子系统和模型窗口的大小和方位


若想查看子系统的内容或对子系统进行再编辑，可以双击子系统模块，则会出现一个显示子系统内容的新窗口。在窗口内除了原始模块外，Simulink自动添加输入模块和输出模块，分别代表子系统的输入端口和输出端口。改变其标签会使子系统的输入输出端口的标签也随之变化。


这里需要注意的是菜单命令Edit→Create Subsystem没有相反的操作命令，即一旦将一组模块封装成子系统，就没有可以直接还原的处理方法了(UNDO除外)。因此，在封装子系统前应将模型保存，作为备份。

(2) 模型验证。

① 模型封装： 我们采用仿真实验的方法在MATLAB的Simulink图形仿真环境下进行模型验证实验，其原理如图326所示。 其中，上半部分为精确模型仿真图，下半部分为简化模型仿真图。



图326模型验证原理图



利用前面介绍的Simulink压缩子系统功能可将验证原理图更加简洁地表示为图327的形式。



图327利用子系统封装后的框图



其中： 






由得到的精确模型和简化模型的状态方程，可得到Fcn、Fcn1、Fcn2和Fcn3的函数形式如下所示： 

Fcn： (0.12*u［1］+0.036*sin(u［3］)*power(u［2］,2)-0.9*sin(u［3］)

*cos(u［3］))/(0.24-0.09*power(cos(u［3］),2))

Fcn1： (0.3*cos(u［3］)*u［1］+0.09*sin(u［3］)*cos(u［3］)*power(u［2］,2)

-6*sin(u［3］))/(0.09*power(cos(u［3］),2)-0.24) 

Fcn2： 0.8*u［1］-6*u［3］

Fcn3： 40*u［3］-2.0*u［1］

② 实验设计： 假定使倒立摆在(θ=0，x=0)初始状态下突加微小冲击力作用，则依据经验知，小车将向前移动，摆杆将倒下。下面利用仿真实验来验证正确数学模型的这一“必要性质”。

③ 编制绘图子程序： 

% Inverted pendulum

% Model test in open loop

% signals recuperation

%将导入xy.mat中的仿真试验数据读出

load xy.mat

t=signals(1,:);%读取时间信号

f=signals(2,:);                      %读取作用力F信号

x=signals(3,:);                      %读取精确模型中的小车位置信号

q=signals(4,:);                      %读取精确模型中的倒摆摆角信号

xx=signals(5,:);                     %读取简化模型中的小车速度信号

qq=signals(6,:);                     %读取简化模型中的倒摆摆角速度信号

% Drawing control and x(t) response signals

%画出在控制力作用下的系统响应曲线

%定义曲线的横纵坐标、标题、坐标范围和曲线的颜色等特征

figure(1)                          %定义第一个图形

hf=line(t,f(:));                      %绘制时间作用力曲线

grid on;%加网格

xlabel('Time(s)')                      %定义横坐标

ylabel('Force (N)')                  %定义纵坐标

axis(［0 1 0 0.12］)                   %定义坐标范围

axet=axes('Position',get(gca,'Position'),...

'XAxisLocation','bottom',...

'YAxisLocation','right','Color','None',...

'XColor','k','YColor','k');   %定义曲线属性

 ht=line(t,x,'color','r','parent',axet);     %绘制时间精确模型小车位置曲线

 ht=line(t,xx,'color','r','parent',axet);    %绘制时间简化模型小车位置曲线

 ylabel('Evolution of the x position')    %定义坐标名称

 axis(［0 1 0 0.1］)                   %定义坐标范围

 title('Response x and x' in meters to a f(t) pulse of 0.1 N')%定义曲线标题名称

 gtext('＼leftarrow f(t)'),gtext('x(t)＼rightarrow')，gtext('＼leftarrow x''(t)')

 % Drawing control and theta(t) response signals

 figure(2)

 hf=line(t,f(:));

 xlabel('Time(s)')

 ylabel('Force(N)')

 axis(［0 1 0 0.12］)

 axet=axes('Position',get(gca,'Position'),...

'XAxisLocation','bottom',...

'YAxisLocation','right','Color','None',...

'XColor','k','YColor','k');

  ht=line(t,q,'color','r','parent',axet);

  ht=line(t,qq,'color','r','parent',axet);

  ylabel('Angle evolution(rad)')

  axis(［0 1-0.3 0］)

  title('Response ＼Theta(t) and ＼theta'(t) in rd to a f(t) step of 0.1 N')

gtext('＼leftarrow f(t)'),gtext('＼theta(t)＼rightarrow'),gtext('＼leftarrow ＼theta''(t)')


④ 仿真实验： 执行该程序之结果如图328所示，从中可见： 在0.1N的冲击力作用下，摆杆倒下(θ由零逐步增大)，小车位移逐渐增加； 这一结果符合前述的实验设计，故可以在一定程度上确认该“一阶倒立摆系统”的数学模型是有效的。同时，由图中我们也可看出： 近似模型在0.8s以前与精确模型非常接近； 因此，也可以认为“近似模型在一定条件下可以表述原系统模型的性质”。




图328模型验证仿真结果




3. 双闭环PID控制器设计 

从图322所示的一阶倒立摆系统动态结构图中不难看出，对象传递函数中含有不稳定的零极点，即该系统为一“自不稳定的非最小相位系统”。

由于一阶倒立摆系统位置伺服控制的核心是“在保证摆杆不倒的条件下，使小车位置可控”(注： 此处本应证明系统的可控性，受篇幅所限，请感兴趣的读者自行证明)； 因此，依据负反馈闭环控制原理，将系统小车位置作为“外环”，而将摆杆摆角作为“内环”，则摆角作为外环内的一个扰动，能够得到闭环系统的有效抑制(实现其直立不倒的自动控制)。

综上所述，设计一阶倒立摆位置伺服控制系统如图329所示。剩下的问题就是如何确定控制器(校正装置)D1(s)/D′1(s)和D2(s)/D′2(s)的结构与参数。



图329一阶倒立摆位置伺服控制系统动态结构图



(1) 内环控制器的设计。

① 控制器结构的选择。

考虑到对象为一非线性的自不稳定系统,故拟采用反馈校正，这是因为其具有如下特点： 

 削弱系统中非线性特性等不希望特性的影响； 

 降低系统对参数变化的敏感性； 

 抑制扰动； 

 减小系统的时间常数。

所以，对系统“内环”采用反馈校正进行控制。如图330为采用反馈校正控制的系统内环框图。




图330反馈控制框图



图中，Ks为伺服电机与减速机构的等效模型(已知Ks=1.6)，反馈控制器D′2(s)可有PD,PI,PID三种形式。那么应该采用什么形式的反馈校正装置(控制器)呢？下面，我们采用绘制各种控制器下的“闭环系统根轨迹”的方法进行分析，以选出一种适合的控制器结构。




图331给出了各种控制器结构下内环系统的根轨迹(其中暂定D2(s)=K为纯比例环节)，从图中不难看出： 采用PD结构的反馈控制器可使系统结构简单，使原来自不稳定的系统稳定。所以，我们选定反馈校正装置的结构为PD 形式的控制器。



图331各种控制器下的内环根轨迹


综上有D′2(s)=KP2+KD2s，同时为了加强对干扰量D(s)的抑制能力，在前向通道上加一个比例环节D2(s)=K。从而有系统内环动态结构如图332所示。






图332系统内环动态结构框图



② 控制器参数的整定。

首先暂定比例环节D2(s)的增益K=-20,又已知Ks=1.6。这样可以求出内环的传递函数为


W2=KKsG2(s)1+KKsG2(s)D′2(s)
=-20×1.6×-2.0s2-401+(-20)×1.6×-2.0s2-40(KP2+KD2s)
=64s2+64KD2s+64KP2-40


由于对系统内环的特性并无特殊的指标要求，对于这一典型的二阶系统采取典型参数整定办法，即以保证内环系统具有“快速跟随性能特性”(使阻尼比ζ=0.7，闭环增益K=1即可)为条件来确定反馈控制器的参数
KP2和KD2，这样就有
  
64KP2-40=6464KD2=2×0.7×64


由上式得

KP2=1.625KD2=0.175
 
系统内环的闭环传递函数为

W2(s)=64s2+11.2s+64


③ 系统内环的动态跟随性能指标。

首先进行理论分析。系统内环的动态跟随性能指标如下： 

固有角频率ωn=64=8

阻尼比ζ=0.7

超调量σ%=e-ζπ1-ζ2×100%=4.6%

调节时间ts≈3ζωn=0.536s(5%允许误差所对应的ts)

下面进行仿真实验。

根据得到的内环系统的闭环传递函数，很容易搭建Simulink仿真模型，如图333所示。




图333搭建的Simulink仿真图



编写绘图子程序如下： 

%将导入simu.mat中的仿真试验数据读出

load simu.mat

t=signals(1,:);

x=signals(2,:);

hf=line(t,x(:));

figure(1)

axis(［0 2 0 1.2］)

grid on

xlabel('Time(s)')

ylabel('The response of the step signal(100%)')

title('Response')

得到的仿真图形如图334所示。从仿真图中可以得知，其响应时间和超调量与理论分析的值相符合。



图334单位阶跃信号作用下的响应曲线



(2) 系统外环控制器设计。

外环系统前向通道的传递函数为

W2(s)G1(s)=64s2+11.2s+64×-0.4s2+10s2
=64(-0.4s2+10)s2(s2+11.2s+64)


可见，系统开环传递函数可视为一个高阶(四阶)且带有不稳定零点的“非最小相位系统”。为了便于设计，需要先对它进行一些必要的简化处理(否则，不便利用经典控制理论与方法对其进行设计)。

① 系统外环模型的降阶［1］。

对于一个高阶系统，当高次项的系数小到一定程度时，其对系统的影响可忽略不计。这样可降低系统的阶次，以使系统得到简化。

首先，对内环等效闭环传递函数进行近似处理。由上可知，系统内环闭环传递函数为

W2(s)=64s2+11.2s+64

若可以将高次项s2忽略，则可以得到近似的一阶传递函数： 
            
W2≈6411.2s+64=10.175s+1

近似条件可以由频率特性导出，即


W(jω)=64(jω)2+11.2(jω)+64
=64(64-ω2)+11.2jω≈6464+11.2jω



所以，近似条件是： ω2c≤6410,即ωc≤2.52。


其次，对象模型G1(s)进行近似处理。

我们知道，G1(s)=-0.4s2+10s2，如果可以将分子中的高次项(-0.4s2)忽略，则环节可近似为二阶环节，即G1(s)≈10s2。

同理，近似条件是： 0.4ω2c≤1010,即ωc≤1.58。

经过以上处理后，系统开环传递函数被简化为

W2(s)G1(s)≈57s2(s+5.7)

近似条件为ωc≤min(2.52,1.58)=1.58。

②  控制器设计［1］。

图335给出了系统外环前向通道上传递函数的等效过程，从最终的简化模型不难看出： 这是一个“Ⅱ型系统”。鉴于“一阶倒立摆位置伺服控制系统”对抗扰性能与跟随性能的要求(对摆杆长度、质量的变化应具有一定的抑制能力，同时可使小车有效定位)，我们可以将外环系统设计成“典型Ⅱ型”的结构形式。同时，系统还应满足前面各环节的近似条件，即系统外环的截止角频率ωc≤1.58。




图335模型简化过程



为了满足以上对系统的设计要求，不难发现所需要加入的调节器D1(s)也应为PD的形式。设加入的调节器为
D1(s)=KP(τs+1)； 同时，为使系统有较好的跟随性能，我们采用单位反馈(D′1(s)=K=1)来构成外环反馈通道，如图336所示。则系统的开环传递函数为

W(s)=W2(s)G1(s)D1(s)=57s2(s+5.7)KP(τs+1)




图336闭环系统结构图



为保证系统剪切角频率ωc≤1.58，不妨取ωc=1.2。对于典型Ⅱ型系统(如图337(a)所示)，其频率特性有如下关系(如图337(b)所示)： 






(1) h=T1T2=5时，为典型Ⅱ型系统最优参数，则有h=τ15.7=5，即τ=0.877，不妨取τ=1。则有系统开环传递函数W(s)=10KP(s+1)s2(0.175s+1)。

(2) K=ω1ωc，则10KP=ω1ωc=1×1.2，即KP=0.12。

至此，图336中的控制器均已求出： D1(s)=0.12(s+1)，D′1(s)=1，D2(s)=-20，D′2(s)=0.175s+1.625。



图337典型Ⅱ型系统频率特性图


综上所述，有图338所示的系统动态结构图。




图338系统动态结构图



4. 仿真实验

综合上述内容，有图339所示的Simulink仿真系统结构图。需要强调的是： 其中的对象模型为精确模型的封装子系统形式。



图339Simulink仿真框图 


系统仿真绘图子程序及仿真结果如下： 

(1) 画图子程序。

% Inverted Pendulum PID

% signals recuperation

%将导入PID.mat中的仿真试验数据读出

load PID.mat

t=signals(1,:);

q=signals(2,:);

x=signals(3,:);

% drawing x(t) and theta(t) response signals

%画小车位置和摆杆角度的响应曲线

figure(1)

hf=line(t,q(:));

grid on            

xlabel('Time(s)')       

ylabel('Angle evolution(rad)')   

axis(［0 10-0.3 1.2］)       

axet=axes('Position',get(gca,'Position'),...

'XAxisLocation','bottom',...

'YAxisLocation','right','Color','None',...

'XColor','k','YColor','k');

 ht=line(t,x,'color','r','parent',axet);

 ylabel('Evolution of the x position(m)')

 axis(［0 10-0.3 1.2］)

 title('＼theta(t) and x(t) Response to a step input')

 gtext('＼leftarrow x(t)'),gtext('＼theta(t)＼uparrow')

其中： 

Fcn1： (0.12*u［1］+0.036*sin(u［3］)*power(u［2］,2)-0.9*sin(u［3］)

*cos(u［3］))/(0.24-0.09*power(cos(u［3］),2))

Fcn2： (0.3*cos(u［3］)*u［1］+0.09*sin(u［3］)*cos(u［3］)*power(u［2］,2)

-6*sin(u［3］))/(0.09*power(cos(u［3］),2)-0.24)


(2) 仿真结果。

图340给出了仿真实验结果。从中可见，双闭环PID控制方案是有效的。



图340系统仿真结果图



为检验控制系统的鲁棒性能，还可以改变倒立摆系统的部分参数来检验系统是否具有一定的鲁棒性。例如，将倒立摆的摆杆质量改为1.1kg，此时的Simulink仿真框图仍为图339，只是将Fcn1和Fcn2作如下修改： 


Fcn1： (0.132*u［1］+0.0436*sin(u［3］)*power(u［2］,2)-1.090*

sin(u［3］)*cos(u［3］))/(0.2772-0.1090*power(cos(u［3］),2))

Fcn2： (0.33*cos(u［3］)*u［1］+0.109*sin(u［3］)*cos(u［3］)*power(u［2］,2)

-6.93*sin(u［3］))/(0.109*power(cos(u［3］),2)-0.2772)



其仿真结果如图341所示。



图341变参数时系统的仿真结果 



从仿真结果可见，控制系统仍能有效地控制并保持倒摆直立，并使小车移动到指定位置，系统控制是有效的。

为了进一步验证控制系统的鲁棒性能，并便于进行比较，我们不妨改变倒立摆的摆杆质量和长度多作几组试验。部分仿真实验的结果见图342和图343。

 

图342摆杆长度不变而摆杆质量变化时系统仿真结果





图343摆杆质量不变而摆杆长度变化时系统的仿真结果



可见，所设计的双闭环PID控制器在系统参数的一定变化范围内能有效地工作，保持摆杆直立并使小车有效定位，控制系统具有一定的鲁棒性。


5. 结论

(1) 本节从理论上证明了所设计的一阶直线倒立摆双闭环PID控制方案是可行的。

(2) 本节的结果在实际应用时(实物仿真)还有如下问题： 

① 微分控制规律易受噪声干扰，具体实现时应充分考虑信号的数据处理问题； 

② 如采用模拟式旋转电位器进行摆角检测，在实际应用中检测精度不佳； 

③ 实际应用中还需考虑初始状态下的起摆过程控制问题。

(3) 一阶直线倒立摆的控制问题是一个非常典型而具有明确物理意义的运动控制系统问题，对其进行深入的分析与应用研究，有助于提高我们分析问题与解决问题的能力。

3．3一阶并联旋转双倒立摆系统建模与模型分析

一阶并联旋转双倒立摆(Parallel Rotating Double Inverted Pendulum，PRDIP)是一种新型的欠驱动基准系统，其系统结构如图344所示。本节所要研究的主要问题包括两个。①能否实现将两个摆杆同时控制到“平衡位置”(垂直不倒)，同时实现旋转臂转角的位置伺服控制(即，旋转臂在平面内任意角度的转动)。对于一类工程实际问题，“实物系统上的实践需首先从理论上证明其可行”。②对于问题①的回答，我们通常根据经验(或直觉)是难以判断出来的； 因此，本节要在对PRDIP系统进行数学建模的基础上，应用现代控制理论的“能控性与能观性”定理来分析系统的“能控性”。

要回答问题①，即实现PRDIP系统的有效控制，可采用现代控制理论的LQR(线性二次型最优控制)等方法，由于它们均是以“系统模型”已知(准确可靠)为前提的，故本节还要针对所建立的模型进行“模型验证”工作，以保障后续的PRDIP控制系统设计工作的有效性［4446］。



图344一阶并联旋转双倒立摆系统构成图


1． PRDIP系统建模

本节应用“拉格朗日方程”方法来建立PRDIP系统模型。对于复杂的动力学系统建模问题，“分析力学”中的拉格朗日方程在机电一体化系统(如，倒立摆系统、两轮平衡车、轮腿机器人等)的建模中扮演着重要的角色。相对于牛顿力学的分析方法，它更加简洁，普适性更好，能够有效地描述各种复杂系统的动力学行为，为控制系统分析和设计提供了一种有效的数学工具。

一阶并联旋转双倒立摆(PRDIP)的系统结构如图345所示。表31给出了PRDIP系统中的基本物理量、含义与单位。  



图345PRDIP系统结构示意图



表31一阶并联旋转双倒立摆系统的主要参数



物理量
含义
单位
m
旋转臂的质量
kg
m1
长摆杆的质量
kg
m2
短摆杆的质量
kg
l1
长摆杆质心到其转轴的距离
m
l2
短摆杆质心到其转轴的距离
m
l
旋转臂长度的一半
m
θ
旋转臂角度
rad
θ1
长摆杆角度
rad
θ2
短摆杆角度
rad
J
旋转臂的转动惯量
kg·m2
J1
长摆杆的转动惯量
kg·m2
J2
短摆杆的转动惯量
kg·m2
g
重力加速度
m/s2
c
旋转臂的摩擦系数
N·m·s/rad
c1
长摆杆的摩擦系数
N·m·s/rad
c2
短摆杆的摩擦系数
N·m·s/rad


在实际的机械系统中，总是存在“齿轮间隙、摩擦、传输延迟”等非线性因素，其将影响实际系统的控制效果。为了便于系统的理论分析，这里仅考虑理想情况下对系统性能影响较大的参数。

在以下数学建模过程中，假设“一阶并联旋转双倒立摆系统”满足以下条件： 

① 仅考虑两个摆杆的转动摩擦，其他形式的摩擦和阻尼均可以认为是微小的参量，可以忽略。

② 将长短摆杆均视为质量分布均匀的刚体细杆，忽略旋转臂与长短摆杆之间连接部分的弹性影响。

③ 驱动电机对于旋转臂的作用是无延时的，且驱动电机总是可以准确无误地输出期望的转矩(或者角加速度)。

④ 旋转臂对长摆杆和短摆杆的作用是无延时的。

⑤ 旋转臂在水平面内做圆周运动，长短摆杆分别在与旋转臂垂直的竖直平面内做圆周运动。

⑥ 旋转臂相对于其旋转轴而言是对称的，即其质心在其旋转轴上。

假设一阶并联旋转双倒立摆系统遵循上述的理想条件，则倒立摆系统的数学建模得到相当大程度的简化，同时也便于理解系统的基本动力学原理，为系统的控制算法设计提供便利。

下面应用“分析力学”的拉格朗日方程，对满足上述条件的一阶并联旋转双倒立摆系统进行数学建模，建模步骤如下。

在表31的基础上，假设旋转臂的质心坐标为(x，y，z)，长摆杆的质心坐标为(x1，y1，z1)，短摆杆的质心坐标为(x2，y2，z2)。假设倒立摆的系统控制输入u为施加在旋转臂上的转矩(或角加速度)，单位是N·m； 倒立摆中的旋转臂、长摆杆、短摆杆角度分别为θ、θ1、θ2。基于图345的坐标系可推导得出旋转臂、长摆杆、短摆杆的质心如下： 

xyz=000(31)
x1y1z1=
lcosθ-l1sinθ1·sinθ
lsinθ+l1sinθ1·cosθ
l1cosθ1
(32)
x2y2z2=
-lcosθ+l2sinθ2·sinθ
-lsinθ-l2sinθ2·cosθ
l2cosθ2

(33)

由摆杆质心的位置可推测出质心的速度

v2i=x·2i+y·2i+z·2i


将θ、θ1、θ2代入，可得旋转臂、长摆杆、短摆杆的质心运动速度为

v2=x·2+y·2+z·2=0(34)
v21=x·21+y·21+z·21=l2θ·2+l21θ·21+2ll1θ·θ·1cosθ1+l21θ·2sin2θ1(35)
v22=x·22+y·22+z·22=l2θ·2+l22θ·22+2ll2θ·θ·2cosθ2+l22θ·2sin2θ2(36)


因为长摆杆和短摆杆是相似的，长摆杆的动能公式也可以应用在短摆杆上，因此这里先考虑长摆杆。

根据上述质心运动速度计算，可得长摆杆的平动动能为

P11=12m1v21(37)


下面考虑长摆杆的转动动能，旋转臂和长摆杆的连接图如图346所示，假设旋转臂的旋转中心为O，旋转臂和长摆杆的连接点为O1。以O为原点、xy轴(为水平面)建立坐标系xyz直角坐标系； 以O1为原点、OO1为x1轴，重力的反方向为z1轴，建立坐标系x1y1z1直角坐标系； 以OO1方向为y2轴，以长摆杆的方向为z轴，建立坐标系x2y2z2。

假设旋转臂的转动角度为θ，长摆杆的转动角度为θ1，O1的线速度(即O1绕O旋转的速度)为v，旋转臂的角速度θ·方向
和v相同。



图346旋转臂和长摆杆的连接图


长摆杆转动动能的计算过程较为复杂，因为角速度是矢量，因此对于长摆杆而言，其转动的角速度不仅仅是长摆杆相对旋转臂转动的角速度θ·1
，旋转臂转动所造成的角速度θ·
使得长摆杆有了“其他”角速度量。对于坐标系x2y2z2来说，长摆杆相对于y2的转动角速度为θ·1
，因为x2、z2、y1处于同一平面，所以相对于x2轴和z2轴的转动角速度是由旋转臂的角速度为θ·
而产生的，相对于x2轴的角速度为θ·

sinθ1
，相对于z2轴的角速度为θ·cosθ1


。

假设长摆杆绕z2的转动惯量为Jz，绕x2和y2的转动惯量为J1，那么长摆杆的转动动能为

P12=12J1θ·21+12J1θ·2sin2θ1+12Jzθ·2cos2θ1(38)


在这里，又因为实际的倒立摆摆杆可近似地视为细长的圆柱体，因此长摆杆以z2为旋转轴时的转动惯量远小于以x2和y2为旋转轴时的转动惯量，即JzJ1，故以下分析忽略Jz。可以认为

P12≈12J1(θ·21+θ·2sin2θ1)(39)


同理，也可以计算得到短摆杆的平动动能和转动动能分别如下

P21=12m2v22(310)
P22≈12J2(θ·22+θ·2sin2θ2)(311)


倒立摆系统的动能由旋转臂、长摆杆和短摆杆的平动动能和转动动能组成，所以倒立摆系统的动能为

K=12mv2+12m1v21+12m2v22+12Jθ·2+12J1(θ·21+θ·2sin2θ1)+12J2(θ·22+θ·2sin2θ2)(312)


假设倒立摆系统的势能为(系统势能最大为零)

P=m1gl1(cosθ1-1)+m2gl2(cosθ2-1)(313)


拉格朗日方程中的拉格朗日函数为系统的动能与势能之差，即

L=K-P=12mv2+12m1v21+12m2v22+12Jθ·2+12J1(θ·21+θ·2sin2θ1)+
12J2(θ·22+θ·2sin2θ2)-m1gl1(cosθ1-1)m2gl2(cosθ2-1)
=12m1(l2θ·2+l21θ·21+2ll1θ·θ·1cosθ1+l21θ·2sin2θ1)+12m2(l2θ·2+l22θ·22+
2ll2θ·θ·2cosθ2+l22θ·2sin2θ2)+12Jθ·2+12J1(θ·21+θ·2sin2θ1)+
12J2(θ·22+θ·2sin2θ2)-m1gl1(cosθ1-1)-m2gl2(cosθ2-1)(314)

根据拉格朗日方程可以得到


ddtLθ·-Lθ=u-cθ·
ddtLθ·1-Lθ1=-c1θ·1
ddtLθ·2-Lθ2=-c2θ·2

(315)


由拉格朗日函数L计算可得

Lθ·=m1l2θ·+m1ll1θ·1cosθ1+m1l21θ·sin2θ1+m2l2θ·+m2ll2θ·2cosθ2+
m2l22θ·sin2θ2+Jθ·+J1θ·sin2θ1+J2θ·sin2θ2(316)
Lθ=0(317)
Lθ·1=m1l21θ·1+m1ll1θ·cosθ1+J1θ·1(318)
Lθ1=-m1ll1θ·θ·1sinθ1+m1l21θ·2sinθ1cosθ1+J1θ·2sinθ1cosθ1+m1gl1sinθ1(319)
Lθ·2=m2l22θ·2+m2ll2θ·cosθ2+J2θ·2 (320)
Lθ2=-m2ll2θ·θ·2sinθ2+m2l22θ·2sinθ2cosθ2+J2θ·22sinθ2cosθ2+m2gl2sinθ2(321)


将以上各式分别代入“拉格朗日方程”算式中，进行归纳整理，并假设2=J2+m2l22
、1=J1+m1l21
、=J+m1l2+m2l2
，即可得到一阶并联旋转双倒立摆的数学模型为


(+1sin2θ1+2sin2θ2)θ¨+m1ll1θ¨1cosθ1+m2ll2θ¨2cosθ2+
1θ·θ·1sin2θ1+2θ·θ·2sin2θ2-m1l1lθ·21sinθ1-m2l2lθ·22sinθ2+cθ·=u
m1ll1θ¨cosθ1+θ¨1-1θ·2sinθ1cosθ1-m1gl1sinθ1+c1θ·1=0
m2ll2θ¨cosθ2+2θ¨2-2θ·2sinθ2cosθ2-m2gl2sinθ2+c2θ·2=0
(322)

为便于控制系统设计，现将上述精确的数学模型在平衡点附近(即，PRDIP系统的双摆均处在垂直朝上的“小角度变化”状态)进行模型的线性化处理(简化的近似处理)； 此时，可以认为sinθi=θi，cosθi=1，可忽略模型中的高次项(如θ·θ·1sin2θ1
等)。则通过整理归纳，我们可以将“线性化”后的PRDIP系统数学模型重新写为“状态方程”表达形式x·=Ax+Bu
，其中

A=
010000
0-I4I5cI6-I2I5m1gl1I6I2I5c1I6-I3I4m2gl2I6I3I4c2I6
000100
0I2I5cI6I7m1gl1I6-I7c1I6I2I3m2gl2I6I2I3c2I6
000001
0I3I4cI6I2I3m1gl1I6I2I3c1I6I8m2gl2I6-I8c2I6
,B=0
I4I5I6
0
-I2I5I6
0
-I3I4I6
x=［θθ·θ1θ·1θ2θ·2］T(323)

I1=J+m1l2+m2l2,I2=m1ll1
I3=m2ll2,I4=m1l21+J1
I5=m2l22+J2,I6=I1I4I5-I22I5-I23I4
I7=I1I5-I23,I8=I1I4-I22(324)


至此，完成了PRDIP系统精确数学模型(微分方程描述)与线性化数学模型(状态空间描述)的建立。下面进行精确数学模型的“模型验证”与基于线性化模型的PRDIP系统“能控性与能观性”分析。

2． PRDIP系统的模型验证

依据上述建模结果(即式(322))，可在MATLAB/Simulink环境中搭建以转矩为输入的一阶并联旋转双倒立摆系统“仿真模型”(相关工作类似2．5．2节与3．2节的“图形化仿真程序”设计方法)。封装后的PRDIP系统Simulink仿真模型如图347所示。

设一阶并联旋转双倒立摆系统的参数如下： 旋转臂的质量为m=1．21kg，长摆杆的质量为m1=0．594kg，短摆杆的质量为m2=0．254kg，设置其他参数J=
0．116kg·m2，J1=0．024kg·m2
，J2=0．0019kg·m2，l取0．4m，l1取0．28m，l2取0．08m，重力加速度g取9．8m/s2，旋转臂、长摆杆、短摆杆的摩擦系数c、c1与c2都取0．0005N·m·s/rad。仿真实验中，设置MATLAB/Simulink仿真为定步长(步长为1×10-5)，求解器选为自动。

下面通过几组实验设计，来验证所建立的数学模型的正确性(设系统的控制输入u为转矩τ)。



图347封装后的PRDIP系统Simulink仿真模型


首先，不考虑系统中存在的摩擦，设定系统输入转矩τ=0，旋转臂的初始角度为θ=0，长摆与短摆的角度初值分别为： θ1=1rad，θ2=πrad。仿真时间设为20s，得到的仿真曲线如图348所示。作为对比，设定长摆与短摆角度初值分别为θ1=-1rad(与前实验1rad关于π对称)，θ2=πrad，其他条件不变，得到的仿真曲线如图349所示。



图348数学模型验证仿真曲线(1)




图349数学模型验证仿真曲线(2)


从图348和图349的仿真结果可知，当不考虑系统摩擦且无输入转矩时，长摆从给定的初始角度处向下运动并开始振荡，整个振荡过程围绕180°位置(即竖直向下位置)进行； 在长摆振荡的影响下，短摆也发生了振荡，整个振荡过程同样围绕180°位置(即竖直向下位置)进行，长摆与短摆的运动过程与实际相符。并且，对比图348和图349的仿真曲线可以发现，由于长摆的初始位置关于竖直方向对称，两次实验的运动过程也是相互对称的，这进一步证明了所建数学模型的正确性。

然后，考虑系统中存在摩擦的情形。设定系统输入转矩τ=0，旋转臂的初始角度为θ=0，长摆与短摆的角度初值分别为θ1=2．5rad，θ2=πrad，仿真时间设为70s，得到的仿真曲线如图350所示。作为对比，分别取θ1=πrad，θ2=2．5rad，其他条件不变，仿真时间设为50s，得到的仿真曲线如图351所示。



图350数学模型验证仿真曲线(3)




图351数学模型验证仿真曲线(4)


从图350的仿真结果可知，当考虑系统中存在的摩擦且输入转矩为零时，长摆从给定的初始位置处开始向下运动并振荡，整个振荡过程围绕180°位置(即竖直向下位置)进行，并且振荡幅度逐渐减小，由于长摆质量重且惯量大，运动过程受摩擦力影响较小，所以长摆用时约70s才停在180°的稳定位置； 在长摆振荡的影响下，短摆也发生了振荡，整个振荡过程同样围绕180°位置(即竖直向下位置)进行，振荡幅度逐渐减小，最终稳定在180°的位置。从图351的仿真结果可知，当考虑系统中存在的摩擦时，由于短摆质量轻且惯量小，运动过程受摩擦力影响较大，所以短摆仅用时约8s就停在了180°的稳定位置； 长摆几乎不受短摆运动的影响，角度波动很小，在±2°范围内。长摆与短摆的运动过程与实际相符，证明了所建数学模型的正确性。

综上所述，在不考虑摩擦与考虑摩擦两种条件下，利用“必要条件法”证明了所建立的PRDIP系统数学模型的正确性。

3．  基于线性化模型的系统能控性与能观性分析

下面采用“控制变量法”仿真分析PRDIP系统参数变化对系统“能控性与能观性”的影响。同时，根据工程经验，我们也发现，长摆与短摆长度之比和系统摩擦系数大小是系统能控性与能观性的主要影响因素，下面一并进行分析。

(1) 长摆与短摆的长度之比和系统能控性与能观性的关系。

在不考虑摩擦的情况下，通过MATLAB软件的rank 语句可计算给出如表32所示的长摆与短摆的长度之比和系统能控性/能观性的关系，表中Q和P分别为系统能控性与能观性矩阵。


表32长摆与短摆的长度之比和系统能控性/能观性的关系


l1/m
l2/m
l1/l2
rank(Q)
rank(P)
0．280
0．020
14
6
6
0．280
0．070
4
6
6
0．280
0．140
2
6
6
0．280
0．280
1
4
6
0．500
0．500
1
4
6
0．100
0．100
1
4
6
0．280
0．560
0．5
6
6
0．280
1．120
0．25
6
6
0．280
2．800
0．1
6
6


由表32中的仿真结果可知，当不考虑摩擦时，无论长摆与短摆长度之比为多少，系统能观性矩阵P均为满秩，即PRDIP系统均是完全能观的； 而当长摆与短摆长度相同时，系统能控性矩阵Q不是满秩，即此时PRDIP系统不是完全能控的； 亦即，只有当长摆与短摆长度不相等时，PRDIP系统才是完全能控的。

(2) 系统摩擦系数大小与系统能控性/能观性的关系。

考虑系统中摆杆存在的摩擦的情况(不考虑旋转臂的摩擦，即C=0)，给出了如表33所示系统摩擦系数大小与系统能控性/能观性的关系。此处，考虑C1=C2。


表33系统摩擦系数大小与系统能控性/能观性的关系


C1=C2/(N·m·s/rad)
rank(Q)
rank(P)
0
6
6
0．0005
6
6
0．0015
6
6
0．005
6
6
0．02
6
6
0．04
6
6
0．05
5
6
0．1
5
6


由表33中的结果可以发现，摩擦系数会影响系统的能控性，当摩擦系数较小时，系统的能控性矩阵才是满秩的，即此时PRDIP系统才可能是完全能控的； 亦即，PRDIP系统中存在“摩擦系数影响系统能控性”问题，故需要在制造实物装置时予以充分考虑。

综上所述，PRDIP系统的参数变化对系统的能观性没有影响，但是对系统的能控性有影响——PRDIP系统是完全能观的； 只有当两个摆杆长度不完全相等且摆杆转动摩擦系数较小时，PRDIP系统才是完全能控的。

4．  结论与思考

通过本节的系统建模与模型分析，我们可以得到以下几点结论。

(1) 本节所建立的PRDIP系统数学模型能够较精确地描述一阶并联旋转双倒立摆系统。其中，线性化之前的模型精度比较高，对于大范围的摆角变化都适用，和实际情况较为接近； 同时，线性化模型是一个“近似模型”，它仅适用于双摆在“平衡点”附近时的“稳摆控制器”设计，亦即便于利用较为成熟的“线性系统理论”进行控制系统设计。

(2) 一阶并联旋转双倒立摆系统是一个自不稳定的“非最小相位系统”，需要通过控制系统设计的方法来实现“闭环控制系统”的稳定； 此处，平衡点处的稳定控制(两摆杆偏离零点不大时)可以基于系统的“线性化模型”进行设计。

(3) 基于PRDIP系统的线性化模型所进行的系统“可控性与能观性”分析可知： 当两根摆杆参数完全一致时候，系统是“不能控”的，也就是在实物系统上不可能实现双摆相同情况下的“稳摆与旋转臂的位置伺服控制”。此外，通过模型分析还发现： 系统旋转部件的“摩擦系数”也会影响系统的可控性，即两摆杆的旋转“摩擦系数”不能太大，否则系统也是不可控的。这一结论为实物实验装置的制造加工提供了有益的参考。

(4) 问题与思考： 针对2．5．2节(龙门吊车运动控制问题)与3．3节(一阶并联旋转双倒立摆系统建模与模型分析)中的拉格朗日方程建模方法，试说明： 为什么所应用的拉格朗日方程数学表达式有所不同？它们之间有何联系？



问题探究
2.5.2节与
3.3节中，
拉格朗日
方程的表
达形式为
何不同？

扫本页的二维码可深入了解。

3.4自平衡式两轮电动车直行与转向复合控制系统设计

现代工业化、航空航天以及人们的生活娱乐不断地为控制技术提出新的问题与挑战，控制理论及其应用就是在探究问题与需求中不断向前发展。随着交通工具向着小型、节能、环保、便携等方向发展，人们开始对微小型“电动车”产生了兴趣。




两轮电动车是继摩托车、汽车之后，人们研制的一种新型代步工具，它仅靠两个轮子来支撑车体，采用电池提供动力，由电动机驱动，采用微控制器、姿态传感器、控制软件及车体机械装置共同协调控制车体的平衡，靠人体重心（或指令）的改变使车辆完成启动、加速、减速、停止等行驶动作。

这个想法刚出现时，人们不禁都要问： 这能行吗？车会不会倒？能否转弯？……

早在2000年，瑞士国家工业电子实验室以Felix Grasser为首的一些研究人员就提出了两轮电动车的设计思想并制作了模型样机（如图352（a）所示）。文献\[2］给出了系统的结构分析、数学建模以及控制方法，开创性地解决了系统的稳定问题，并考虑到干扰情况下的鲁棒性问题。但是，其系统采用的是固定闭环极点配置的线性状态反馈控制方法，使得控制效果并不理想，而且没有给出车体在转弯、上下坡等方面问题的深入研究。

2001年7月，自平衡式两轮电动车产品——Segway HT（human transporter）首先在美国问世，由迪恩·卡门设计（如图352（b）所示）。Segway HT的核心技术是“动态稳定控制技术”，它的工作原理类似于人身体的平衡控制方式（人类身体由内耳、眼睛、大脑、肌肉等来控制平衡），产品的动态稳定控制应用了电晶体陀螺仪、倾角传感器、软件与硬件电路板（微处理器）、电机来工作。



图352自平衡式两轮电动车系统


Segway HT产品的详细情况，可浏览网站http://humantransporter.com/。


为拓展自平衡式两轮电动车的应用领域，适应社会需要，Segway公司于2009年初联手美国通用公司推出了一款坐式Segway概念车——PUMA（如图352（c）所示）。它在外形与时速上更贴近普通机动车，并可两人同时乘坐，周身布置的多个避障传感器，增加了它的安全性能，提高了其实用性。

在国内，中国科学技术大学也研制了此类电动车，采用动力学理论对自平衡式两轮电动车进行力学分析，建立了系统的数学模型，采用状态反馈等控制理论实现了其稳定运行［3］。

2010年，通用汽车在上海首发了以PUMA为基础、融合电气化和车联网两大技术的双座ENV电动联网概念车（Electric NetworkedVehicle）（如图352（d）所示）。该车时速可达40km/h，充电一次可行驶40km，重量为400kg，车身体积1.5m×1.435m×1.64m，自身携带3个GPS可实现精确定位。车身通过电机拖动以实现在平台上的前后滑动，从而改变车身重心位置。ENV是通用汽车对未来城市个人交通的最新解决方案，可使未来城市交通实现零油耗、零排放、零堵塞和零事故。

综上可见，自平衡式两轮电动车（以下简称“自平衡两轮车”）运动控制的主要问题是“直行与转向复合控制问题”，本节将基于拉格朗日方程法建立自平衡两轮车的数学模型，对模型进行验证与分析，推导出直行与转向状态下系统的等效动态模型； 同时，根据系统左右轮输入转矩的对称耦合特性，设计了自平衡两轮车直行与转向复合控制系统； 最后，基于MATLAB/Simulink仿真工具，对自平衡两轮车“直行与转向复合控制方案”进行仿真实验验证，以证明该方案的有效性与可行性。

1. 系统建模

为便于对自平衡两轮车的运动控制方案进行仿真实验分析，以及实际产品的开发，需要建立其有效精确的数学模型； 同时，在不影响系统性能的情况下，需要对系统中难以处理的部分进行适当的假设，以便于控制系统的分析与设计。

我们假设： 车身与驾驶者等效为一根刚性直杆； 车身处于平衡直立状态时，质心位于两轮中心正上方； 运行时忽略车轮与地面间的滑动摩擦。

在建立自平衡两轮车数学模型中，所涉及物理变量及其所对应的物理意义，如表34所示。


表34变量名及其意义



符号
说明
符号
说明
θt
车身倾角
Ml/Mr
左/右轮电机转矩

θl/θr
左/右轮转角
R
车轮半径
L
车身质心到轮轴的距离
D
车轮间距
fl/fr
左/右轮阻尼力
Jt
车身绕轮轴的转动惯量
φ
车身转角
Jz
车身绕Z轴的转动惯量
mb
车身质量
Jw
车轮绕轮轴的转动惯量
mw
车轮质量
Jφ
车轮绕Z轴的转动惯量



为便于分析与计算，建立自平衡两轮车的车体坐标系Oxb,yb,zb与地坐标系Oxe,ye,ze［4］。在车体坐标系Oxb,yb,zb中，令zb轴过自平衡两轮车车身质心并垂直于轮系轴线，由原点Ob指向质心； xb轴为轮系轴线延长线，由原点指向右轮圆心； yb轴垂直于xb轴与zb轴，由原点指向自平衡两轮车正前方。在地坐标系Oxe,ye,ze中，ze轴指向重力加速度g的负方向，xe轴指向正东方向，ye轴指向正北方向。令自平衡两轮车处于初始状态时，两坐标系原点重合。


当车身前倾角度为θt时，将yb轴投影到xeOye平面上，设为y′b轴，则此时y′b轴与ye轴重合。再令自平衡两轮车绕ze偏转φ角度，过mb质心作垂线垂直于xeOye面，交y′b轴于A点，此时即自平衡两轮车运行的一般状态，其物理模型可以简化为图353所示。下面，我们将应用拉格朗日方程法对自平衡两轮车系统进行数学建模。

自平衡两轮车系统共有三个自由度： 车身绕轮轴xb轴的前倾后仰、车身绕ze轴的转动以及车轮在xeOye面内沿y′b轴的前进与后退。而车身的转动是由左右轮的转速差引起的，故选取车身倾角θt、左轮转角θl与右轮转角θr作为广义变量。

（1） 速度与动能计算。

为应用拉格朗日方程法进行建模，首先需求解系统的总动能，车体总动能可以分解为车身与车轮的平动动能和转动动能。

对于车身速度，车身的前倾后仰、旋转以及前行后退都会引起车身速度的变化，下面分三种情况来进行速度求解。

车身前倾后仰时，考虑车身前倾情况，其速度可以分解为图354所示。其中，v1=θ·tL为车身前倾引起的速度大小，则其三个分量的大小可以表示为

vz1=－θ·tLsinθt（325）
vx1=θ·tLcosθtsinφ（326）
vy1=θ·tLcosθtcosφ（327）



图353自平衡两轮车简化物理模型




图354车身前倾速度分解示意图



车身旋转时，取图353的俯视图，其速度可以分解为图355所示。其中，线段OA=Lsinθt，故v2=φ·Lsinθt为车身旋转引起的速度大小，则此速度的两个分量大小可以表示为
vx2=－φ·Lsinθtcosφ（328）
vy2=φ·Lsinθtsinφ（329）

车身前行后退时，考虑其前行情况，取图353的俯视图，其速度可以分解为图356所示。其中v3=vl+vr/2为车身前进引起的速度大小，则其两个分量可以表示为
vx3=12vl+vrsinφ（330）
vy3=12vl+vrcosφ（331）
其中
vl=θ·lR
vr=θ·rR（332）


图355车身旋转速度分解示意图




图356车身前进速度分解示意图



注意： 车身绕Z轴的转动惯量Jz是个不定值，它随摆角θ变化。由于要控制车身的平衡，即θ≈θr（θr为角度给定值，随路面情况变化），因而Jzθ变化范围很小，近似看作常值，即Jzθ≈Jzθr≈Jz。

因此，车身平动速度为
vx=vx1+vx2+vx3
vy=vy1+vy2+vy3
vz=vz1（333）
车身平动动能为
T1=12mbv2x+v2y+v2z（334）
车身转动动能为
T2=12Jtθ·2t+12Jzφ·2（335）
车轮平动动能为
T3=12mwv2l+v2r（336）
车轮转动动能为
T4=12Jwθ·2l+θ·2r+2×12Jφφ·2（337）
其中，Jw=12mwR2，Jφ=mwD22。



图357转弯轨迹分解示意图


在自平衡两轮车数学模型中，为减少变量数量，需利用已有的几何关系，对模型进行简化处理，进而消去中间变量。这里，车身转向角φ可通过两个车轮的转角经过变换得到，取自平衡两轮车俯视图如图357所示，将车身的转向的轨迹MP近似分解为车身绕右轮旋转的轨迹MN与车身直线运动的轨迹NP，即存在如下近似的几何关系
MP≈MN+NP（338）
代入各变量，可以认为
θlR=φD+θrR（339）
整理可得
φ=θl－θrRD（340）
故系统总动能为
T=T1+T2+T3+T4
=12mbθ·2tL2+θ·l－θ·r2L2R2sin2θtD2+θ·l+θ·r2R24+θ·tθ·l+θ·rLRcosθt+
12Jtθ·2t+JzR2θ·l－θ·r22D2+34mwR2θ·2l+θ·2r+14mwR2θ·l－θ·r2（341）

（2） 应用拉格朗日方程建模［47］

对自平衡两轮车进行受力分析,可得系统在三个广义坐标方向上的广义力分别为
Qθt=mbgLsinθt－Ml－Mr（342）
Qθl=Ml－flR（343）
Qθr=Mr－frR（344）
考虑左右轮受到的粘滞阻尼力，定义fl=klθ·l，fr=krθ·r，其中，kl与kr为粘滞阻尼系数。

注意： 根据自平衡式两轮电动车的驾驶工艺，车身倾角θt越大，两轮转速θ·i也越快，反之亦然，可以近似认为θt∝θ·i。当两轮车恒速运行时，两轮车车体满足转矩关系fl+frR=mbgLsinθt≈mbgLθt，因此有fi∝θ·i，即fl=klθ·l，fr=krθ·r，粘滞阻尼系数kl与kr随路面状况不同而改变。

通过以上计算，求得了系统的总动能与广义力，将式（341）、式（342）、式（343）以及式（344）代入拉格朗日方程即可得到系统的精确数学模型，故系统在θt方向上有
mbL2+Jtθ¨t+12mbcosθtLRθ¨l+θ¨r－mbL2R2sinθtcosθtθ·l－θ·r2D2=mbgLsinθt－Ml－Mr（345）
在θl方向上有
12mbRLθ¨tcosθt+mbR2L2sin2θtD2+14mbR2+JzR2D2+2mwR2θ¨l+
－mbR2L2sin2θtD2+14mbR2－JzR2D2－12mwR2θ¨r－12mbRLθ·2tsinθt+
2mbR2L2θ·tθ·l－θ·rsinθtcosθtD2=Ml－klθ·lR（346）
在θr方向上有
12mbRLθ¨tcosθt+mbR2L2sin2θtD2+14mbR2+JzR2D2+2mwR2θ¨r+
－mbR2L2sin2θtD2+14mbR2－JzR2D2－12mwR2θ¨l－12mbRLθ·2tsinθt－
2mbR2L2θ·tθ·l－θ·rsinθtcosθtD2=Mr－krθ·rR（347）
式（345）、式（346）和式（347）即为自平衡两轮车系统的精确数学模型。

（3） 模型线性化。

为了便于应用已有理论对系统进行设计，通常要对复杂的系统模型进行简化。在这里，需要对上述精确数学模型进行线性化处理。考虑到自平衡两轮车在实际运行过程中，θt≤10°，故sinθt≈θt，cosθt≈1，同时忽略一些高次项，如sin2θt≈0，θ·2t≈0，θ·tsinθt≈0，θ·l－θ·r2sinθt≈0。此外，系统受到的阻尼力与电机提供的转矩相比也可以忽略不计，则将精确模型线性化后的模型变为
mbL2+Jtθ¨t+12mbLRθ¨l+θ¨r=mbgLθt－Ml－Mr（348）
12mbRLθ¨t+14mbR2+JzR2D2+2mwR2θ¨l+14mbR2－JzR2D2－12mwR2θ¨r=Ml （349）
12mbRLθ¨t+14mbR2+JzR2D2+2mwR2θ¨r+14mbR2－JzR2D2－12mwR2θ¨l=Mr （350）
（4） 模型验证。

建立了自平衡两轮车的数学模型后，需要对其进行仿真验证，判断其是否符合实际，选取实际物理参数值如表35所示。


表35数学模型参数取值



参数
取值
参数
取值
mb
70kg
R
0.2m
mw
6kg
D
0.5m
L
0.7m
Jt
36.62kg·m2
Jz
1.3kg·m2
kl
9.5N·s/rad
g
9.8kg·m/s2
kr
9.5N·s/rad


这里采用“必要条件法”来检验模型。所谓“必要条件法”，是指所进行的模型验证实验结果是依据经验可以预知的，其正确的结果是“正确的模型”所应具备的“必要性质”。根据实际情况，将上述参数代入精确数学模型可得

5.488sin2θt+1.388θ¨l+0.372－5.488sin2θtθ¨r+
4.9θ¨tcosθt－θ·2tsinθt+10.976θ·tθ·l－θ·rsinθtcosθt=Ml－1.9θ·l（351）
5.488sin2θt+1.388θ¨r+0.372－5.488sin2θtθ¨l+
4.9θ¨tcosθt－θ·2tsinθt－10.976θ·tθ·l－θ·rsinθtcosθt=Mr－1.9θ·r（352）
70.92θ¨t+4.9cosθtθ¨l+θ¨r－5.488θ·l－θ·r2sinθtcosθt
=480.2sinθt－Ml－Mr（353）

为验证所建立的数学模型是否正确，在MATLAB/Simulink环境下进行仿真以确定是否符合实际。在Simulink中搭建的模型如图358所示。




图358自平衡两轮车数学模型仿真结构图


系统模型封装如图359所示。其中两个输入量Ml、Mr分别代表左右轮电机的输入转矩。Fcn1、Fcn2和Fcn3中的函数为由式（351）、式（352）和式（353）经过整理得出的θ¨t、θ¨l、θ¨r的表达式。


Fcn1： 

(-((-5.488*sin(u(5))*sin(u(5))+0.372)*u(4)+10.976*(u(1)-u(3))*u(6)*sin(u(5))*cos(u(5))+4.9*(u(7)*cos(u(5))-u(6)*u(6)*sin(u(5))))+u(2)-1.9*u(1))/(5.488*sin(u(5))*sin(u(5))+1.388)

Fcn2： 

(-((-5.488*sin(u(5))*sin(u(5))+0.372)*u(4)+10.976*(u(1)-u(3))*u(6)*sin(u(5))*cos(u(5))+4.9*(u(7)*cos(u(5))-u(6)*u(6)*sin(u(5))))+u(2)-1.9*u(1))/(5.488*sin(u(5))*sin(u(5))+1.388)

Fcn3： 

(4.9*cos(u(1))*(u(7)+u(6))-5.488*(u(4)-u(5))*(u(4)-u(5))*sin(u(1))*cos(u(1))-480.2*sin(u(1))+u(2)+u(3))/(-70.92)


六个输出量中，Bs代表车身倾斜角速度，Ba代表车身倾角，Ls与Rs分别代表左右轮转动角速度，La与Ra分别代表左右轮转角，以上各量单位为弧度制。



图359自平衡两轮车数学模型封装图



下面利用两组仿真实验来对模型进行验证。

实验一，通过修改仿真模型中积分器的初始值，给定自平衡两轮车的车身初始倾角为10°，即θt=0.35rad，取输入转矩Ml=Mr=0N·m。考虑实际情况中，在不受两轮电机控制的情况下，车身会朝向初始倾角方向倾倒，车轮朝反方向滚动，最终趋于静止。仿真实验结果如图360所示。



图360模型验证实验一仿真结果与实际情况


从实验一仿真结果可以看出，车身倾角经过20s左右最终平衡在180°附近，即车身竖直向下。这是由于在仿真中没有对车身倾斜角度的约束，即没有实际情况中地面的阻挡，使得车身可以指向正下方，即偏转了180°。因此，实验一中，仿真结果与实际情况相符。


实验二，在初始时刻给定自平衡两轮车左轮电机15N·m初始转矩，右轮电机10N·m初始转矩，车身保持直立，即θt=0°，Ml=15N·m，Mr=10N·m。考虑实际情况中，车身会朝车轮运行方向的反方向倾倒，在两轮电机的驱动作用下，车轮持续向前滚动，其中，左轮转角会大于右轮转角，车体实际上进行左转弯运行。实验二仿真结果如图361所示。



图361模型验证实验二仿真结果与实际情况


从实验二仿真结果可以看出，车身倾角朝车轮运动方向的反方向倾倒，约20s后最终平衡于-180°，即方向竖直向下。而在两轮电机的驱动作用下，左轮转角大于右轮转角。因此，实验二中，仿真结果与实际情况相符。

综合以上两个仿真实验，仿真结果与“实际经验”相一致； 因此，可以认为仿真实验在一定程度上验证了前文所建立模型的正确性。

（5） 直行与转向运行系统等效动态结构图。

为便于进行直行与转向复合控制系统的设计，基于建模给出如下直行与转向两种情况的系统等效动态结构图。

当自平衡两轮车直行时，其左右两个车轮可等效为一个车轮的运动，此时，令θl=θr=θw，θ¨l=θ¨r=θ¨w，Ml=Mr=Mw，代入线性化模型式（348）、式（349）与式（350）中，整理可得
mbL2+Jtθ¨t+mbRLθ¨w=－2Mw+mbgLθt
12mbRLθ¨t+12mbR2+32mwR2θ¨w=Mw（354）
将上述模型验证所取实际参数代入式（352）可得
70.92θ¨t+9.8θ¨w=－2Mw+480.2θt

4.9θ¨t+1.76θw=Mw（355）
对式（355）进行拉普拉斯变换可得
70.92s2－480.2θts+9.8s2θws=－2Mws
4.9s2θts+1.76s2θws=Mws（356）
故系统传递函数模型为
G1s=θtsMws=－0.17s2－10.99G2s=θwsθts=－6.13s2+36.73s2（357）
系统等效动态结构图如图362所示。



图362直行系统等效模型


将图362与本书3.2节中图322所介绍的一阶直线倒立摆模型相比，可看出它们具有相同的结构。因此，自平衡两轮车直行控制策略可参考一阶直线倒立摆的控制策略。

当两轮以轮轴中心为圆心原地旋转时，此时两轮轨迹为圆形，并认为车身处于相对静止状态，为方便分析，令θt=0°。此时，θl=－θr，θ¨l=－θ¨r，Ml=－Mr，代入线性化模型式（348）、式（349）与式（350）中，整理可得
2JzR2D2+4mwR2θ¨l=Ml（358）
带入实际参数，进行拉普拉斯变换得到系统传递函数模型为
Gls=θlsMls=11.376s2（359）
Grs=θrsMrs=11.376s2（360）


图363转向运行时的系统
等效模型

因此，转向运行时的系统等效模型如图363所示。


2. 自平衡两轮车运动控制系统的控制器设计

（1） 自平衡两轮车系统能控性与能观性分析。

自平衡两轮车是一个自不稳定系统，要实现对其有效的控制，应首先验证其“能控性与能观性”； 因此，我们先要对系统的能观性/能控性进行分析。

自平衡两轮车的数学模型经过线性化，带入实际参数可得
70.92θ¨t+4.9θ¨l+4.9θ¨r=－Ml－Mr+480.2θt
4.9θ¨t+1.388θ¨l+0.372θ¨r=Ml
4.9θ¨t+0.372θ¨l+1.388θ¨r=Mr（361）
为了方便将上式表示成状态方程，选取θt、θ·t、θ·l和θ·r作为状态变量，将式（361）表示成矩阵形式为
1000070.924.94.904.91.3880.37204.90.3721.388θ·tθ¨tθ¨lθ¨r=0100480.200000000000θtθ·tθ·lθ·r+00－1－11001MlMr
（362）
化成标准状态方程形式为
θ·tθ¨tθ¨lθ¨r=010011.005000－30.638000－30.638000θtθ·tθ·lθ·r+00－0.087－0.0871.0180.0330.0331.018MlMr（363）
因此，令x=θtθ·tθ·lθ·r，y=θ·lθ·r，u=MlMr，故由上式可得到三个矩阵分别为A=

010011.005000－30.638000－30.638000，B=00－0.087－0.0871.0180.0330.0331.018，C=00100001

即可将式（361）表示成x·=Ax+Bu
y=Cx的标准方程形式。

能控性定理告诉我们： 如果一个系统是能控的，就可以对这个系统进行控制器设计，以实现有效地控制。

系统的能控矩阵可表示为
M=BAB…An-1B（364）
根据MATLAB语言中ctrb函数的计算，可得到上述自平衡两轮车系统的能控矩阵为
M=00－0.087－0.08700－0.9574－0.9574－0.087－0.08700－0.9574－0.9574001.0180.033002.61592.6159000.0331.018002.61592.615900
（365）

通过计算得到rankM=n=4，依据能控性定理知： 自平衡两轮车系统是可控的，车身倾角和车轮转速均可通过设计适当的控制器加以控制。

能观性定理告诉我们： 如果一个系统是能观的，就可以对这个系统进行状态观测器设计，以实现对状态变量的有效测量(直接或间接的)。

系统的能观矩阵可表示为
N=CCA︙CAn-1（366）
根据MATLAB中obsv函数的计算，可以得到上述自平衡两轮车系统的能观矩阵为
N=00100001－30.638000－30.6380000－30.638000－30.63800－337.1712000－337.1712000（367）

通过计算得到rankN=n=4，说明自平衡两轮车系统是可观测的； 即，系统中如有无法直接测量的状态变量，可以通过设计适合的状态观测器来进行“软测量”。

（2） 直行状态下系统的双闭环PID控制器设计。

由前推导，我们已经得到了自平衡两轮车系统直行时的等效动态结构图如图362所示。由于车轮角速度与车轮转角存在ωs=sθws的关系，则可将动态结构图表示成如图364所示。从图364中可以看出，系统传递函数中含有不稳定的零极点，即系统是自不稳定的“非最小相位系统”。



图364直行系统动态结构图


对于自平衡两轮车的直行运动，控制系统的控制目标有二： 车身倾角（直立不倒），车轮速度（快慢可调），而前者又是后者的基础； 因此，自平衡两轮车直行运动控制采用“双闭环控制方案”，将倾角环设计为内环，速度环设计为外环。综上，自平衡两轮车直行运动时的双闭环控制系统动态结构图如图365所示。



图365自平衡两轮车直行运动时的双闭环控制系统动态结构图


鉴于自平衡两轮车在直行时的模型与一阶直线倒立摆系统的线性化模型在结构形式上是一致的，基于“类比原理”，我们采用本书3.2节中的“一阶直线倒立摆双闭环PID控制方案”设计方法，设计自平衡两轮车直行时的双闭环PID控制器参数。其中，倾角环控制器参数为
D2s=－423.66（368）
D′2s=1.15+0.17s（369）
为使转速环取得最佳抗扰性能，我们令D1s为PI控制器，使系统为典型Ⅱ型系统，可有外环控制器设计参数如下： 
D1s=0.019+0.0075s（370）
D′1s=1（371）
根据以上设计所得参数，对自平衡两轮车直行双闭环PID控制策略的有效性在MATLAB/Simulink下进行仿真实验验证。搭建系统仿真模型如图366所示。



图366直行控制器仿真模型图


编写绘图子程序如下： 

clf

load shiyan3_.mat

t = signals(1,:);

Ls = signals(2,:);

Rs = signals(3,:);

Ba = signals(4,:);

Lst = line(t,Ls(:));

Rst = line(t,Rs(:));

grid on;

xlabel('t/s');

ylabel('车轮转速/(m/s)');

axis(［0 30 -2 2］);

axet=axes('Position',get(gca,'Position'),'XAxisLocation','bottom','YAxisLocation','right','Color','None','XColor','k','YColor','k');

Bat = line(t,Ba(:),'color','k','parent',axet);

ylabel('车身倾角/deg');

axis(［0 30 0 8］);

gtext('＼leftarrow 车身倾角'),gtext('＼leftarrow 车轮转速')

设车身初始倾角为0°，在初始时刻给定车轮转速为1m/s。仿真结果如图367所示。



图367直行控制器仿真结果


从仿真图形可以看出，车轮转速经过小幅度超调，在10s左右稳定在给定的1m/s，此时自平衡两轮车做匀速直线运动，车身为了保持平衡，倾角经小幅负调，最终稳定在2.2°左右，这与实际相符。

综上可见，仿真实验验证了自平衡两轮车直行时双闭环PID控制策略的有效性。

（3） 直行与转向复合控制系统设计。

前文已经推导出自平衡两轮车线性化后的数学模型式（348）、式（349）与式（350），为便于下面的分析，现令
a=mbL2+Jt（372）
b=12mbLR（373）
c=14mbR2+JzR2D2+2mwR2（374）
d=14mbR2－JzR2D2－12mwR2（375）
e=mbgL（376）
则以上三式变为
aθ¨t+bθ¨l+bθ¨r=eθt－Ml－Mr（377）
bθ¨t+cθ¨l+dθ¨r=Ml（378）
bθ¨t+dθ¨l+cθ¨r=Mr（379）
联立消去θ¨l与θ¨r可得
ac+d－2b2θ¨t=ec+dθt－b+c+dMl+Mr（380）
进行拉普拉斯变换可得
θts=－b+c+dac+d－2b2s2－ec+dMls+Mrs（381）
此外，对式（378）与式（379）联立，整理可得
θ¨l=1c2－d2cMl－dMr－bc+dθ¨t（382）
θ¨r=1c2－d2cMr－dMl－bc+dθ¨t（383）
进行拉普拉斯变换，并代入ωls=sθls，ωrs=sθrs得到
ωls=1c2－d2scMls－dMrs－bsc+dθts（384）
ωrs=1c2－d2scMrs－dMls－bsc+dθts（385）
由式（381）、式（384）与式（385）可以得到自平衡两轮车系统动态结构图，如图368所示。其中，Gs=－b+c+dac+d－2b2s2－ec+d。




图368自平衡两轮车系统动态结构图



由图368可以看出，系统两驱动轮（输出转矩）之间存在着“耦合关系”（相互影响），当左/右轮输出（控制设定）大小相等方向相同的转矩时，自平衡两轮车沿直线行驶； 当左/右轮输出（控制设定）大小相等方向相反的转矩时，自平衡两轮车原地转向运行。可以推断： 自平衡两轮车在平面内的任何运动形式均可以由直线运行（两轮输出转矩相等）与原地转向运行（两轮输出转矩相反）合成产生； 因此，在实现自平衡两轮车的平面运动控制中，我们采取“分别独立设计直行与转向控制器”的方法，以避开系统所存在的复杂的“耦合问题”（理论上也可以采用“解耦控制”方法），实现系统的有效控制［5］。

图369给出了基于图368数学模型的仿真结果，其中： 忽略系统阻尼，初始时刻车身倾角为0°，给定左轮转矩为2N·m，右轮转矩为0N·m。从仿真实验结果不难看出，左轮给定初始转矩后，右轮也随之一同转动，自平衡两轮车左右轮间存在耦合关系。



图369耦合存在性仿真实验结果


下面我们来设计自平衡两轮车的直行与转向复合控制系统。

首先，对于直行控制部分，图365已给出“双闭环PID控制方案”，仿真实验证明其可行，因此这里我们依然采用之； 依据图368的系统耦合模型，可建立直行控制系统动态结构如图370所示（上半部分），其中为求得车体运行的平均速度，在速度环反馈中进行了算术平均处理； 直行控制部分的系统输入可以看作自平衡两轮车的“油门”。


其次，对于转向控制部分，我们设： Mws为两轮平均转矩，ΔMws为两轮转矩差，即Mws=Mls+Mrs/2，ΔMws=Mls－Mrs； 则依据“转向等价转矩差”思想，有如图370所示（下半部分）的转向控制系统动态结构，其中Δωs为两轮车转向运行控制器输入，输入给定量为两轮转速差ωl－ωr，可以看作两轮车的“方向盘”。



图370自平衡两轮车直行与转向复合控制系统动态结构图


在图370转向控制器D3s设计中，基于前述图363的“转向运行状态系统传函”，有如图371所示的转向控制系统结构。



图371自平衡两轮车转向运行控制系统结构图


在转向环设计中，为保证转向控制在阶跃/斜坡给定下实现无差跟踪，系统应具有Ⅱ型结构； 因此，令D3s为PI控制器，采取单位负反馈，即D′3s=1。

设
D3s=Kp3+Ki3s（386）
则转向环闭环传递函数为
W3s=Kp3s+Ki31.376s2+Kp3s+Ki3（387）
由式（387）可见，转向环闭环传递函数为“具有零点的二阶系统”，相对于常规的二阶系统其具有更好的跟踪快速性［6］。

将W3s表成典型形式为
W3s=Ki31.376Kp3Ki3s+1s2+Kp31.376s+Ki31.376（388）
由“二阶系统最佳参数”理论［6］，令阻尼比ζd=0.707，调节时间ts=0.5s(误差带±5%)，可列出如下方程
2ζdωn=Kp31.376
ω2n=Ki31.376
ts=3ζdωn=0.5（389）
其中，ωn为系统固有角频率。

由上可解得转向环控制器D3s参数为
D3s=16.51+99.1s（390）


根据上述直行与转向复合控制系统方案，设计自平衡式两轮电动车载人驾驶/操控系统如图372所示。



图372自平衡式两轮电动车载人驾驶/操控系统结构图


其中，驾驶者通过人工视觉获取自平衡两轮车的车速信息与转弯角度信息，通过油门与方向盘给定自平衡两轮车的车速大小与行驶方向； 同时，车身自动保持平衡控制，符合人们常规的驾驶习惯。由上述设计可见，基于系统建模的直行与转向复合控制系统设计方案，使得自平衡两轮车的操控与人们的驾驶习惯相一致，为自平衡两轮车的安全驾驶提供了有效保障。


3. 仿真验证

根据以上计算所得控制器参数，对控制系统在MATLAB/Simulink下进行仿真实验验证，搭建系统模型如图373所示。



图373自平衡两轮车控制系统仿真模型图



仿真实验参数给定如表36所示。


表36仿真实验参数给定值



时刻/s
0
15
25
40
50
油门给定/(m/s)
0.8
0.8
0.8
0.8
0.8
方向盘给定/(m/s)
0
0.4
0
-0.4
0




图374自平衡两轮车理论行驶轨迹


由表36可知，自平衡两轮车行驶的平均速度给定始终为0.8m/s，在第15s时方向盘向右旋转，两轮车向右转向，25s时方向盘回正，两轮车直线行驶，40s时方向盘向左旋转，两轮车向左转向，50s时两轮车再次直线行驶，仿真实验中，左右两轮理论行驶轨迹如图374所示，仿真结果如图375、图376所示。绘图子程序如下： 

clf

load shiyan3_.mat

t = signals(1,:);

Ls = signals(2,:);

Rs = signals(3,:);

Ba = signals(4,:);

Lst = line(t,Ls(:));

Rst = line(t,Rs(:));

grid on;

xlabel('t/s');ylabel('车轮转速/(m/s)');

axet=axes('Position',get(gca,'Position'),'XAxisLocation','bottom','YAxisLocation','right','Color','None','XColor','k','YColor','k');

Bat = line(t,Ba(:),'color','k','parent',axet);

ylabel('车身倾角/（°）');

gtext('＼＼leftarrow 车身倾角'),gtext('＼＼leftarrow 车轮转速')，gtext('＼＼leftarrow 两轮行驶距离')



图375仿真实验速度曲线




图376仿真实验行驶距离与倾角曲线


由图375与图376可见，两轮车初始时刻直线运行，两轮行驶距离相等，15s时两轮车向右转向运行，左轮行驶距离开始大于右轮，25s时直线行驶，40s时两轮车向左转向运行，两次转向行驶时间与速度大小均相同； 因此，从图376中可以看出50s后两轮行驶距离再次相等，与实际相符； 同时，从倾角仿真曲线可以看出，车身倾角在初始时刻最大为4.5°左右，行驶时保持在1.8°左右。

综上，整个行驶过程中车身保持平衡，两轮速度与给定相符，仿真结果验证了直行与转向复合控制系统设计方案的有效性。

4. 小结

本节基于自平衡式两轮电动车系统的数学模型，给出了车体直行与转向复合控制方案，通过仿真实验验证了系统设计方案的有效性。

另外，为使两轮车能适应不同身高/体重的人驾驶，不但要求控制系统能有效地实现车体的平衡、直行与转向运行控制，还要保证整车系统在上下坡、加减速、恶劣路况等条件下，仍然具有良好的操控性（即鲁棒性），这也是自平衡式两轮车得以实用化的关键所在。受篇幅所限，本书不便深入，感兴趣的读者可以在本节内容的基础上展开深入研究。


3.5龙门吊车重物防摆的滑模变结构控制方案

1. 问题提出


在前面一节中，我们应用鲁棒PID控制方案实现了固定绳长条件下，龙门吊车的重物防摆与系统的定位控制。然而，吊车在实际运行中绳长常常需要变化，我们把吊车实际的工作过程抽象为图377所示情况，其“搬运—定位—安放”过程为由A点提升至B点(或由B点下放至A点)。




图377吊车重物定位过程示意图



显然，前面一节所述方法在解决这类问题时就无能为力了。因此，需要研究一种新的控制方法，以实现变绳长条件下的重物防摆与系统的定位控制。在这里将采用“滑模变结构控制方案”，实现变绳长条件下龙门吊车重物防摆与定位控制。



同时，为检验滑模变结构控制方法的有效性，将使用MATLAB/Simulink软件工具，对所设计的控制系统进行仿真实验研究。由于滑模变结构控制是一种“不连续的控制方法”，系统在进行Simulink仿真时将会面临一些特殊的问题。因此，本节中还将探讨一类具有“不连续控制率”的控制系统的Simulink仿真方法。


2. 滑模变结构控制［79］

变结构控制理论诞生于20世纪50年代末。作为一种非线性控制理论，与其他控制器相比，它具有控制规律简单，对系统的数学模型精确性要求不高，可以有效地调和动、静态之间的矛盾以及具有强鲁棒性等优点。近年来，已被广泛应用于处理一些复杂的非线性、时变、多变量耦合及不确定系统的控制中，如伺服电机驱动、机器手控制以及飞行器控制系统的设计等。下面对滑模变结构控制的一些基本概念作以简要阐述。

(1) 滑动模态。

滑模变结构控制是变结构控制系统(可简称为VSS)的一种控制策略。这种控制策略与常规控制的根本区别在于“控制的不连续性”，它是一种“系统结构随时变化的开关特性”，该控制特性可以迫使系统在一定条件下沿规定的状态轨迹作小幅度、高频率的运动，称为滑动模态或“滑模”运动。

这种滑动模态是可以设计的，且与系统的参数及扰动无关。这样，处于滑模运动的系统就可以具有很好的鲁棒性。

(2) 切换函数与切换面。

在滑模变结构控制中，需要通过开关的切换，改变系统在状态空间中的切换面s(x)=0两边的结构。开关切换的法则称为控制策略，它保证系统具有可滑动的模态。此时，分别把s=s(x)及s(x)=0叫做切换函数及切换面。

(3) 滑模变结构控制的数学描述。

滑模变结构控制可表述成如下形式。

设有一非线性控制系统： 

x·=f(x,u,t)x∈Rn，u∈Rm，t∈R

确定其切换函数向量： 

s(x),s∈Rn

具有的维数一般情况下等于控制的维数。

寻求变结构控制率： 

ui(x)=u+i(x)当si(x)＞0
u-i(x)当si(x)＜0




图378滑模变结构控制的基本原理

这里的变结构体现在u+(x)≠u-(x)，使得



① 满足“到达条件”： 切换面si(x)=0以外的相轨迹将于有限时间内到达切换面； 

② 切换面是滑动模态区，且滑动运动是渐进稳定的，动态品质应良好。满足以上条件的控制称为“滑模变结构控制”，其基本原理可用图378表示。


(4) 滑模变结构控制系统的设计方法。

滑模变结构控制系统的设计步骤可概括如下。

① 选择滑模面参数，构成希望的滑动模态； 

② 求取不连续控制u±，保证在切换平面s=0上的每一点存在滑动模态，这一平面就被认为是滑动面； 

③ 控制必须使状态进入滑动面。

滑动模态三要素(即存在、稳定、进入)是靠切换面与控制两者来保证的，一旦切换面选定，则问题转入控制的求取。求取标量变结构控制，要从滑动模态的存在条件出发，即从ss·<0这个关系式出发。按这个关系式所求得的控制，往往是不等式。在选取时，可充分考虑滑动模态的进入条件。对于滑动模态的稳定性问题，由于在选择切换面时已经考虑过，因此不需要再进行讨论。


3. 系统建模及模型验证

在2.5.2节中已经讨论了吊车系统的模型建立和模型验证问题，这里不再赘述，下面将着重讨论吊车系统滑模变结构控制问题。


4. 基于滑模变结构控制的吊车重物防摆控制策略

这里我们主要研究变绳长条件下的吊车重物防摆控制，在2.5.2节的式(244)中，给出了变绳长吊车系统的数学模型(如下式)： 

(m0+m)x¨1-mx¨2sinx3-mx2x¨3cosx3-2mx·2x·3cosx3+
mx2x·23sinx3+Dx·1=f1mx¨2-mx¨1sinx3-mx2x·23-mgcosx3=f2mx22x¨3+2mx2x·2x·3-mx¨1x2cosx3+mgx2sinx3+ηx·3=0

忽略空气阻尼的影响η=0，可将上面的方程组转换成如下形式，以便于滑模变结构控制器的设计： 


x¨1=1m0(f1-Dx·1+f2sinx3)
x¨2=gcosx3+x2x·23+f1-Dx·1m0sinx3+m0+msin2x3m0mf2
x¨3=-gx2sinx3-2x·2x2x·3+cosx3x2m0(f1-Dx·1+f2sinx3)



其中，m0和m分别是小车和重物的质量，x1、x2和x3是系统的状态，分别代表小车的位置、绳长和摆角，f1和f2是系统的输入，分别是小车受到的水平方向的拖动力和绳长方向受到的拉力，D为轮轨摩擦系数，g为重力加速度(取9.8m/s2)。

设x1=x1，x2=x2，x3=x3，x4=x·1，x5=x·2，x6=x·3，则可得系统状态方程为

x·1=x4
x·2=x5
x·3=x6
x·4=(f1-Dx4+f2sinx3)/m0
x·5=gcosx3+x2(x6)2+sinx3(f1-Dx4)/m0+f2(m0+msin2x3)/m0m
x·6=-gsinx3/x2-2x5x6/x2+(f1-Dx4+f2sinx3)cosx3/m0x2


在此数学模型的基础上，将设计一个多输入向量滑模变结构控制器，使吊车既可实现精确定位，又能实现绳长达到期望值，同时保证摆角最小。

(1) 定义两个滑模面： 

s1=x·1+α1(x1-P)+α3x3+α4x·3=x4+α1(x1-P)+α3x3+α4x6s2=x·2+α2(x2-L)=x5+α2(x2-L)

其中，α1,α2均是正实数。

(2) 采用等速趋近率： 

dsdt=-ρsgn(s)

则可保证ss·<0，满足广义滑模条件： 


s·1=-W1sgn（s1）s·2=-W2sgn（s2）
W1>0,W2>0


经过计算整理，可以得到

x·4=-α4x·6-α1x4-α3x6-W1sgn(s1)

x·5=-α2x5-W2sgn(s2)
x·6=-2x5x6-gsinx3-α1x4cosx3-α3x6cosx3-W1cosx3sgn(s1)x2+α4cosx3


(3) 求控制力： 

吊车系统的电机水平拉力f1为

f1=m0x·4+Dx4-(sinx3)f2

电动机提升力f2为

f2=-m(sinx3)x·4+mx·5-mx2(x6)2-mgcosx3

由此可知，水平拉力f1和电动机提升力f2均为系统的全状态反馈，需要控制系统的各个状态向量，则由系统状态方程式以及控制力f1和f2即可构成变绳长吊车防摆定位滑模变结构控制器。


5. 仿真实验

(1) 具有“不连续控制率”的系统仿真方法研究［10，11］。

上面所求出的控制力的表达式中存在不连续性函数sgn(s)，这给系统仿真带来了特殊问题。

为了分析说明的简便，下面以固定绳长防摆的滑模变结构控制为例，介绍“不连续控制率”的系统仿真方法。这时，控制方式为只调节内环摆角，令x1=θ，x2=θ·，l=1m，a=x¨，这里x、l和θ分别代表小车的位置、绳长和摆角，则在线性区附近，系统方程可简化为如下二阶状态方程： 


x·1=x2x·2=-9.8x1+a


按照前面介绍的滑模变结构控制系统的设计方法，可以得到具有防摆功能的二维滑模变结构控制器为a=9.8x1-c1x2-ρsgn(c1x1+x2)。

这里我们全部采用Simulink工具箱中的模块搭建系统的仿真模型，如图379所示。




图379二维滑模变结构控制器仿真模型



取c1=1，ρ=1，x1(0)=1，x2(0)=0，仿真时间T=10s，采用Variablestep的ODE45算法。仿真结果如图380所示。




图380系统时域仿真结果



由图380可以看出，当系统仿真到1s左右时，速度非常慢，仿真停滞不前。

结果分析： 这是由于系统存在不连续模块sgn，当系统于1s左右到达滑模面(s=0)时，sgn模块向系统发出过零通知。而当采用变步长求解器时，Simulink能够检测到过零。所以，当Simulink检测到过零时，便自动缩小步长，可是在下一仿真步长里，系统继续过零。因为滑模面在1s处不能正常归零，所以sgn模块就反复过零，同时一直向求解器发出过零通知。求解器便相应地一直不停地缩小步长，这样由于仿真步长太小，系统便在不连续处形成了过多的点，超出了系统可用的内存和资源，使得系统进展缓慢，仿真停滞不前。

基于以上原因，我们提出了以下4种解决策略： 

① 取消Zero crossing detection功能； 

② 采用fixedstep求解器； 

③ 采用不能够产生过零通知的Fcn函数模块； 

④ 柔化sgn(s)函数，使其连续化。


以上所提出的4种解决策略中，采用单独任何一种或几种都可行。这里我们采用第三种方法，使用Fcn函数模块来解决仿真停滞问题，图381为按照这种方法搭建的仿真模型。



图381采用Fcn函数模块后的系统仿真模型



图中，Fcn的表达式为： -9.8*u［2］+u［3］； 

Fcn1的表达式为： 9.8*u［2］-u［3］*u［1］-u［4］*sgn(u［3］*u［2］+u［1］)。

取c1=1，ρ=1，系统仿真结果如图382所示。




图382使用Fcn函数模块后的系统时域仿真结果



由于Fcn函数模块不支持过零检测，所以系统在不连续的情况下仍然能迅速完成仿真。除此之外，采用fixedstep求解器和变步长下置Zero crossing detection为off的仿真方法，其仿真结果亦如图382所示。其实这三种方法的本质相同，都是取消了过零检测。由上面的仿真结果可以看出： 

① 取消系统的过零检测功能之后，仿真速度快； 

② 但是，由于仍然存在不连续模块sgn，所以加速度存在较大的抖振。

对于第四种方法： 柔化sgn(s)函数策略，既可以解决仿真停滞问题，又可消除加速度的抖振，是解决这一问题的最佳方法，在下文中我们将对此详细阐述。

在解决了仿真停滞问题的基础上，下面讨论更符合吊车实际运行情况的变绳长防摆定位控制系统的仿真。


(2) 变绳长吊车防摆滑模变结构控制系统的仿真。

根据前面设计的控制器表达式，在MATLAB/Simulink环境下建立变绳长吊车防摆的滑模变结构控制系统仿真模型，如图383所示，这里采用Fcn函数模块来解决上面提出的仿真停滞问题。




图383变绳长情况下吊车防摆定位滑模变结构控制系统模型



取x1~x6，f1，f2，m0，m，D分别作为输入u［1］~u［11］。则根据吊车系统的数学模型，可得Fcn0，Fcn1和Fcn2的函数形式分别为

Fcn0：  (u［7］-u［11］*u［4］+u［8］*sin(u［3］))/u［9］

Fcn1： 9.8*cos(u［3］)+u［2］*u［6］*u［6］+(u［7］-u［11］*u［4］)*sin(u［3］)/

u［9］+u［8］*(u［9］+u［10］*sin(u［3］)*sin(u［3］))/(u［9］*u［10］)

Fcn2： -9.8*sin(u［3］)/u［2］-2*u［5］*u［6］/u［2］+(u［7］+u［8］*

sin(u［3］)-u［11］*u［4］)*cos(u［3］)/(u［9］*u［2］)

取x1~x6，α1~α4，W1，W2，P，l，m0，m，D分别作为输入u［1］~u［17］。则根据控制力f1和f2的表达式，即可得Fcn3和Fcn4的函数形式为 

Fcn3： u［17］*u［4］+u［15］*(-u［7］*u［4］-u［9］*u［6］-u［11］*
sgn(u［4］+u［7］*(u［1］-u［13］)+u［9］*u［3］+u［10］*u［6］)
-u［10］*((-2*u［5］*u［6］-9.8*sin(u［3］)-u［7］*u［4］*
cos(u［3］)-u［9］*u［6］*cos(u［3］)-u［11］*cos(u［3］)*sgn(u［4］+u［7］*(u［1］-u［13］)+u［9］*u［3］+u［10］*u［6］))/(u［2］+u［10］*cos(u［3］))))-sin(u［3］)*(-u［16］*u［2］*u［6］*u［6］-9.8*
u［16］*cos(u［3］)-u［16］*sin(u［3］)*(-u［7］*u［4］-u［9］*u［6］
-u［11］*sgn(u［4］+u［7］*(u［1］-u［13］)+u［9］*u［3］+u［10］*u［6］)-u［10］*((-2*u［5］*u［6］-9.8*sin(u［3］)-u［7］*u［4］*
cos(u［3］)-u［9］*u［6］*cos(u［3］)-u［11］*cos(u［3］)*sgn(u［4］+u［7］*(u［1］-u［13］)+u［9］*u［3］+u［10］*u［6］))/(u［2］+u［10］*cos(u［3］))))+u［16］*(-u［8］*u［5］-u［12］*sgn(u［5］+u［8］*(u［2］-u［14］))))

Fcn4： -u［16］*u［2］*u［6］*u［6］-9.8*u［16］*cos(u［3］)-u［16］*
sin(u［3］)*(-u［7］*u［4］-u［9］*u［6］-u［11］*sgn(u［4］+u［7］*(u［1］-u［13］)+u［9］*u［3］+u［10］*u［6］)-u［10］*((-2*u［5］*u［6］-9.8*sin(u［3］)-u［7］*u［4］*cos(u［3］)-u［9］*u［6］*
cos(u［3］)-u［11］*cos(u［3］)*sgn(u［4］+u［7］*(u［1］-u［13］)+u［9］*u［3］+u［10］*u［6］))/(u［2］+u［10］*cos(u［3］))))+u［16］*(-u［8］*u［5］-u［12］*sgn(u［5］+u［8］*(u［2］-u［14］)))


为使仿真模型更加简洁清晰，可以使用Simulink封装子系统功能对吊车模型进行封装，得到模型如图384所示。




图384封装后的滑模变结构控制系统模型



图384模型中的crane system模块为封装后的吊车系统模型。打开crane system模块，可得吊车系统模型如图385所示。其输入为电机水平拉力f1和电机提升力f2，输出为吊车的位置、绳长和摆角。




图385吊车系统模型




对于滑模变结构控制器，取α1=17，α2=15，α3=1，α4=-0.45，W1=20，W2=5，P=3m，l=1m，m0=50kg，m=10kg，D=0.1； 对于初始状态，取x2(0)=2m，x1(0)=x3(0)=x4(0)=x5(0)=x6(0)=0。仿真结果如图386所示，这里分别绘制了小车位置、绳长变化、摆角大小、水平拉力和电机提升力运行轨迹曲线。从图中可以看出，小车位置、绳长在2.5s内迅速归于期望值，与此同时，摆角也归于零点，并稳定下来。但是，水平拉力f1和提升力f2出现了高频抖振现象，使得系统难以工程实现。因此，希望通过改进算法来降低系统抖振。




图386滑模变结构控制系统的仿真曲线



(3) 准滑模伪变结构控制系统的仿真研究。

① 抖振的产生。

当具体实现理想变结构系统时，理想的开关特性u(x)=u*(x)sgn［s(x)］是不可能实现的。时间延迟和空间滞后等因素将使得滑动模态表现为抖动形式，在光滑的滑动上叠加了自振。这种现象称为抖振。

抖振问题是变结构控制广泛应用的突出障碍，是影响变结构技术发展的重要原因。这是因为： 有的系统元件不能承受高频切换； 有的系统性能上不允许存在抖振； 抖振的存在还可能激发系统未建模部分的强迫振荡。于是在实践中，人们尝试采用各种具有准滑动模态的控制系统。所谓“准滑动模态”(或近似滑动模态、伪滑动模态)，是指系统的运动轨迹被限制在理想滑动模态的某一Δ邻域内的模态。

② 抖振的抑制。

抖振发生的本质原因是开关的切换动作造成控制的不连续性。因此，对一个现实的滑模变结构控制系统，抖振必定存在。我们可以努力去削弱抖振的幅度，使它减少到工程允许的范围内，而无法完全消除它，消除了抖振，也就消除了滑模变结构控制系统的抗干扰能力。由于抖振问题是滑模变结构控制的突出障碍，许多学者提出了消抖措施，其中有代表性的有前面介绍的柔化sgn(s)函数法、边界层法和趋近率法等。这里采用“柔化sgn(s)函数法”，即用一连续函数U(s)

U(s)=ss+δ≈ss=sgn(s)


来代替不连续函数sgn(s)，如图387所示。这种降低抖振的方法，也被称为“继电函数连续化的准滑模伪变结构控制”。由于将继电函数连续化，所以系统不再存在结构变化，故称为“伪变结构”。




图387继电函数及继电函数连续化曲线




③ 准滑模伪变结构控制。

采用“柔化sgn(s)函数法”的准滑模伪变结构控制算法来降低抖振后，系统不再存在抖振，且各状态变量效果良好。但是，系统控制力仍然较大，远远超出了实验系统所能承受的范围(大于1000N)，究其原因，乃滑行速度过大所致。所以，进一步将原滑行速度降低至W1=1.5，W2=0.5。系统的仿真结果如图388所示(取δ=0.1)，可见使用柔化sgn(s)函数法既解决了仿真停滞问题，又消除了抖振，是解决不连续控制率系统仿真问题的最佳方法。



图388消振后的实验结果



(4) 准滑模伪变结构控制系统的鲁棒性研究。

从理论角度讲，滑动模态是降维光滑运动，且渐近趋向原点，而且系统的滑模运动与控制对象的参数变化、系统的外部扰动及内部的摄动无关，因此滑模变结构控制系统的“鲁棒性”要比一般常规的连续控制系统强得多。

下面，研究继电函数连续化的变绳长吊车防摆定位准滑模伪变结构控制器的鲁棒性能。如图389~图392所示，它们分别为在不同的载荷质量、不同的吊车自重、不同的摩擦系数以及外部扰动情况下的系统响应曲线。



图389载荷质量不同时的实验结果




图390吊车自重不同时的实验结果





图391摩擦系数不同时的实验结果





图392正弦函数扰动前后的实验结果



通过上面的仿真实验结果可以看出:

① 系统在内部参数(载荷质量、吊车自重、摩擦系数)变化时，位移、绳长以及摆角都能保持原运动轨迹不变； 

② 系统在外部扰动下，位移、绳长以及摆角都能在有限的时间内归于期望值，且无较大的波动。

6. 结论

(1) 本节设计了基于滑模变结构控制的变绳长吊车防摆控制器，较为系统地介绍了滑模变结构控制理论，总结归纳了滑模变结构控制系统的设计方法。它可概括为： 首先，选择适当的滑模面参数，构成希望的滑动模态； 其次，求取不连续控制u±，保证在切换平面s=0上的每一点存在滑动模态，并使系统状态进入滑动平面。


(2) 本节探讨了具有“不连续控制率”的系统仿真方法，解决了仿真的停滞问题。在MATLAB/Simulink仿真时，由于控制力中存在不连续性函数sgn(s)，如果直接使用Simulink工具箱中的sgn模块，当发生过零检测时仿真将会停滞。为此我们提出了四种解决这一问题的方案，其中使用继电函数连续化的准滑模伪变结构控制，既解决了仿真停滞问题，又消除了抖振现象，是解决不连续控制率系统仿真问题的最佳方法。

这也提示我们，在使用MATLAB/Simulink进行系统仿真时，对于工具箱中的一些模块，常常需要根据实际情况进行合理选择和设计，以使得仿真更符合系统实际情况，且在保持仿真精度的前提下提高仿真速度。

(3) 仿真结果表明，所设计的滑模变结构控制系统具有很强的鲁棒性。通过在不同载荷质量、不同吊车自重、不同摩擦系数以及外部扰动情况下的系统鲁棒性仿真实验，可以得出滑模变结构控制在内部参数变化和外部扰动下，具有比其他一般控制方法更强的鲁棒性，便于实际应用。

3．6基于可视化模型的PRDIP运动控制系统设计
1． 可视化模型

对于一般的机电一体化系统建模，可先基于动力学原理建立系统的精确数学模型，以备用于控制系统的设计； 但是，对于一类复杂的机电一体化系统，例如一阶并联旋转双倒立摆系统(PRDIP)这样的多连杆、强欠驱动、强耦合的非线性机械系统，基于动力学原理(如，分析力学的拉格朗日方程)的数学建模的过程是相当复杂的，而且如果对模型的某个局部环节分析稍不合理，将会导致系统动力学建模结果的不佳，进而影响基于系统数学模型的控制系统设计与性能分析结果的不准确，甚至出现严重错误。

随着计算数学理论与人工智能技术的不断进步，基于各种先进算法的自动建模技术与“工具软件产品”应运而生； 在这个领域比较有代表性就是MathWorks公司开发的MATLAB/Simscape工具箱，连同相关专业模块等辅助产品，MATLAB/Simscape已成为现今应用广泛的数字仿真工具软件。在针对复杂机电一体化系统的建模工作中，MATLAB/Simscape可有效避免出现上述数学建模过程的繁杂与不确定性，对于像一阶并联旋转双倒立摆系统(PRDIP)一类的复杂工程系统，可以采用多领域物理建模的方法，建模过程和装配实际硬件系统一样，将描述硬件元件的软件模块按照一定规律逐一连接起来建立仿真模型，而系统所基于的数学模型会在组装过程中“自动建模”出来。

(1) 可视化模型与MATLAB/Simscape(Multibody)工具箱。

Simulink： Simulink是集成在MATLAB中的一个可用于动态系统建模和仿真的环境。用户不仅可以利用其中的模块库和求解器等进行建模和仿真，还可将MATLAB算法融入Simulink模型中进行联合仿真，仿真结果还可导出到MATLAB中进行分析。

Simscape： Simscape是Simulink的一个模块，它是一种面向对象的多领域物理仿真工具，其中的模块集除了基础模块库(Foundation Library)，还包括集成度更高的电气与电子模块集(Electrical)、动力传动系统模块集(Driveline)、多体机械系统模块集(Multibody)等。用户可以使用该工具箱提供的这些模块集在Simulink环境中迅速创建物理系统的模型。

Simscape Multibody： Simscape Multibody是Simscape的一个模块，称为多体机械系统模块集(Multibody)，其中包括刚体子模块组(Body Elements)、关节模块组(Joint)、Force and Torque(力与转矩模块组)、Frames and Transforms(坐标系与变换模块组)和Constraints(约束模块组)等。用户可以使用该模块集提供的模块组对多体机械系统进行建模，而系统所基于的数学模型会在组装过程中自动建立起来，同时用户还可以在自动生成的3D动画中查看系统动态。

关于应用Simscape Multibody工具箱更详细的建模方法与应用案例，可查阅MATLAB中相关的Simscape工具箱使用方法以及网上的一些可视化建模实例视频教程。

(2) 基于MATLAB/Simscape工具箱的PRDIP系统建模。

采用Multibody模块集中的Cylindrical Solid(固体圆柱)作为PRDIP系统的旋转臂、摆杆、底座和支柱，并根据实际需要设置尺寸大小和材料。采用Multibody模块集中的Revolute Joint(旋转关节)和Rigid Transform(刚体变换模块)，作为旋转臂两端与摆杆的连接部件。在Revolute Joint(旋转关节)中可以设置摆杆初始角度和转动摩擦系数。在Rigid Transform(刚体变换模块)中设置合适的参考系以限制摆杆的旋转平面。

根据以上分析，在MATLAB/Simulink仿真环境下，搭建了如图393所示的仿真模型。



图393基于MATLAB/Simscape工具箱的PRDIP系统可视化模型仿真图


图393中，Solid1为PRDIP底座，Solid2为支撑柱，Solid3为旋转臂，Solid4为长摆，Solid5为短摆，Solid4的内部如图394所示，可设置其尺寸与材料等参数； RT1、RT2、RT3、RT4和RT5均为刚体变换环节，其主要作用为调整该环节前后两项的相对连接位置，使其符合设计要求，RT1内部如图395所示； RJ1、RJ2和RJ3均为转动关节，其中B端为输入端，F端为输出端，摆杆初始角度和摩擦等相关参数可在其内部进行设置，RJ2内部如图396所示。



图394Solid4内部图





图395刚体变换环节RT1内部图




图396旋转关节RJ2内部图





图397PRDIP系统可视化模型图

运行后得到了如图397所示的一阶并联旋转双倒立摆系统可视化模型。可视化模型的建模过程更加简便，建模思路便于理解，得到的模型更加直观而有效，实现了PRDIP系统仿真模型动态变化过程的“可视化”。结合图393形成的图397的可视化模型，同样可以和3．3节建立的数学模型一样进行仿真模型的封装，对于图393的可视化Simulink模型，形成的封装形式如图397所示(对应图347的Simulink封装形式)。这样，我们可在对两种建模方式设置同样的参数后，进行两种方式建模的准确性对比验证； 即，在后续的控制设计中直接调用图397的可视化模型进行控制算法的仿真实验分析。

2．基于可视化模型的模型分析

作为商业化的工具软件，MATLAB/Simscape工具箱具有较高的可信度； 同时，作为检验数学模型的准确性与有效性这一“模型验证”工作，在“必要条件法”的模型验证基础上，所建立的数学模型与现代仿真工具得到的可视化模型的“对比分析”，也是一种可行的“模型验证”方法。因此，此处将建立的可视化模型与3．6节中建立的数学模型的仿真结果进行对比，以此证明所建立的PRDIP系统可视化模型与数学模型的正确性与一致性。

值得注意的是，在对比仿真验证过程中要保证可视化模型和数学模型中的参数设置完全一致，一般可根据实际装置系统的参数建立可视化模型，基于可视化模型中的参数设置进一步获得数学模型中的一些参数，在相同参数下进行仿真实验和对比分析。

首先，考虑系统中存在的摩擦，设定系统输入转矩τ=0，旋转臂的初始角度为θ=0，长摆与短摆的角度初值分别为θ1=2．5rad，θ2=πrad，仿真时间设为70s，得到的仿真曲线如图398所示。



图398考虑摩擦时可视化模型仿真结果


将图350所示的数学模型仿真结果与图398所示的可视化模型仿真结果进行作差处理，得到了如图399所示的摆杆角度差值曲线。



图399数学模型与可视化模型仿真结果差值曲线


然后，同样考虑系统中存在的摩擦，设定系统输入转矩τ=0，旋转臂的初始角度为θ=0，长摆与短摆的角度初值分别为θ1=πrad，θ2=2．5rad，仿真时间设为50s，得到的仿真曲线如图3100所示。



图3100可视化模型验证仿真实验摆杆角度仿真曲线


将图351所示的数学模型仿真结果与图3100所示的可视化模型仿真结果进行作差处理，得到了如图3101所示的摆杆角度差值曲线。



图3101数学模型与可视化模型仿真结果差值曲线


通过对比图350和图398、图351和图3100，以及观察图399与图3101可知，所建立的PRDIP系统数学模型与可视化模型，在同一初始条件得到的仿真结果基本上是一致的(整体误差小于1°，在可接受范围内)，证明了本节所建立的PRDIP系统的数学模型和可视化模型是正确的、一致的、可用的。

3．   PRDIP系统的运动控制方案设计与仿真实验

一阶并联旋转双倒立摆(PRDIP)系统全域运动控制过程示意图如图3102所示，这里的“全域运动控制”是指： 起摆控制(摆角可在360°范围内变化)、稳摆控制(摆角小于或等于±1°）与旋转臂的位置伺服控制(旋转臂可在360°范围内任意位置角度的转动)。考虑到摆杆可在其运动平面四个象限内摆动，本节将PRDIP系统的全域控制定义为： 摆杆起摆控制+摆杆稳摆控制+旋转臂位置伺服控制。因此，PRDIP系统全域控制的目标为： 在双摆自然下垂的初始状态下，通过自动控制实现双摆自动摆起到直立不倒的平滑切换控制，同时实现旋转臂位置的精确伺服控制。



图3102PRDIP系统全域运动控制过程示意图


(1) 基于能量反馈控制策略的倒立摆起摆控制方案设计。

PRDIP系统的起摆控制可分为“串行起摆”(先长后短)与“并行起摆”(长短同时)两种方式，本节采用“串行起摆”控制方案。

起摆控制方案主要包括BangBang控制方案、能量反馈控制方案，以及逆系统轨迹控制方案等，考虑到摆杆临界状态的控制效果和控制方案的实现难易，本节选择了能量反馈控制控制方案。在实际起摆控制中，根据已有经验，一般先控制长摆(摆杆1)起摆，起摆成功后稳定摆杆1的同时，再起摆短摆(摆杆2)。先考虑单独摆起摆杆1时，摆杆2对于摆杆1而言是一个干扰项，那么当然希望这个干扰项越小越好； 因此我们假设，摆杆2的摆长质量相对于摆杆1的摆长质量而言较小，这样就可以完成起摆模态的第一步，即摆杆1的自起摆。进而当摆杆1起摆直至进入稳定控制区域范围内时，单独切换摆杆1 的控制模态为稳摆模态，此时摆杆2的运动对于摆杆1 而言仍是一个较小的干扰项。下一步对摆杆2进行起摆控制，由于此时摆杆1已经在稳摆模态中，所需要的电机转矩并不大，而摆杆2相较摆杆1的摆长和质量都较小，其所需要的起摆转矩也应较小，所以此时，选取一个合适的起摆转矩，便可实现既能保证摆杆1不倒，又可以实现摆杆2 的摆起，最终当摆杆2摆起至稳定控制区域范围内时，切换控制器至双摆稳摆模态。

从图3102可以直观地看出，PRDIP系统的起摆控制就是通过自动控制将摆杆从竖直向下位置摆到接近于竖直向上位置。这个过程中，摆杆的位置发生了变化，从较低的位置摆到了较高的位置，摆杆的能量是不断增加的，说明只要给摆杆输入足够的能量，它就能够实现自动摆起。下面以长摆的起摆控制器设计过程为例说明起摆控制器的设计思路。

首先，规定长摆杆处于竖直向上位置时且静止时候的其能量E1ref为0，即为参考能量。考虑长摆一般状态时，采用3．3节规定的PRDIP系统变量符号，则可以得到长摆具有的能量包括动能和势能E1：  

E1=J1θ·21/2+m1gl1(cosθ1-1)(391)


为了避免烦琐的求解微分方程的过程，并且能够判断一阶并联旋转双倒立摆系统的稳定性，可采用设计李雅普诺夫函数的方法来对一阶并联旋转双倒立摆起摆问题进行分析。由于需要对起摆过程中的摆杆能量进行控制，可选取如下李雅普诺夫函数： 

V1=12(E1-E1ref)2+12λ1Jθ·2,λ1>0(392)

式中，除去传统的能量二次方部分，还引入了旋转臂角速度的平方项，这是因为起摆过程中需要对摆杆的角速度和旋转臂的转速进行限制，因为起摆过程的最终目标还涉及与稳摆模态控制器的切换。对于简易的倒立摆，例如一阶直线倒立摆，其稳摆控制器比较成熟，鲁棒性抗干扰性能都较强，其对于起摆过程中的角速度限制要求不高，而此处的控制对象更为复杂，需要对起摆过程的状态变量进行限制，主要体现在旋转臂的转速和摆杆的角速度上。但因为一阶并联旋转双倒立摆系统是一个具有强欠驱动特性的系统，无论是旋转臂的转速还是摆杆的角速度，最后都取决于系统的输入量，即电机输出转矩。换言之，电机转矩直接与旋转臂的转速相关，因此在李雅普诺夫函数的选取上，加入旋转臂角速度的平方项，即对于角速度的限制项，以免因为旋转臂的超调过大引起系统的失稳。

对式(392)的李雅普诺夫函数进行求导可得

dV1dt=E1(J1θ·1θ¨1-m1gl1sinθ1θ·1)+λ1Jθ·θ¨(393)


通过对摆杆1的受力力矩进行平衡分析可以得到

J1θ¨1=Jθ¨cosθ1+m1gl1sinθ1(394)

唯一的控制力矩作用在旋转臂上形成加速度，此处记为

τout=Jθ¨(395)

将式(394)和式(395)代入式(393)并整理可以得到

dV1dt=τout(θ·1E1cosθ1+λ1θ·)(396)


由于能量函数为摆杆1的总能量和旋转臂动能的平方和，另外又规定摆杆1位于竖直向上的稳定状态时能量为0，所以摆杆1在起摆的过程中，能量函数应逐渐减小至0。为保证能量损失尽可能小，因此能量函数随时间的变化应该是单调递减的，即式(396)应始终小于0，即能量函数的导数小于0，保证了摆杆1起摆过程中的能量是逐渐趋近于0，摆杆1的状态是由自然下垂趋近于竖直向上的。为保证式(396)始终小于0，可取如下转矩表达式：  

τout=-k1(θ·1E1cosθ1+λ1θ·),k1>0(397)

从而使得李雅普诺夫函数的导数始终为负，即

dV1dt=-k1(θ·1E1cosθ1+λ1θ·)2≤0,k1>0(398)

从而确保了能量函数是单调递减的，即摆杆1的运动状态是趋近于竖直向上的。

进一步为了保证起摆的快速性，可以令系统的控制输入始终为某一最大值，仅改变其方向。引入符号函数sgn，那么可令控制量为

τout=-k1sgn(θ·1E1cosθ1+λ1θ·),k1>0(399)


短摆起摆时，可类似长摆的推导过程，得到系统所受力矩可取为如下形式： 

τout=-k2sgn(θ·2E2cosθ2+λ2θ·),k2>0,λ2>0(3100)

式中，类似的设定短摆位于竖直向上位置时所具有的能量E2ref为0，E2为短摆所具有的能量，其表达式如下：  

E2=J2θ·22/2+m2gl2(cosθ2-1)(3101)


至此，基于能量反馈控制策略的起摆控制器设计完成。下面将通过仿真实验来验证所设计的起摆控制器的有效性。

在考虑PRDIP系统中存在摩擦的情况下，在MATLAB/Simulink仿真环境下，搭建了如图3103所示的基于PRDIP可视化模型的仿真模型，以此验证基于能量反馈控制策略的起摆控制器对长摆的控制效果。在PRDIP可视化模型的旋转关节(Revolute Joint)中，设定长摆的初始角度θ1=πrad=180°，短摆的初始角度θ2=πrad=180°，即二者初始时均处于自然下垂状态，单独对长摆进行起摆控制，得到了如图3104所示的仿真结果。



图3103长摆起摆控制仿真模型图




图3104长摆起摆控制仿真结果


由图3104的仿真结果可知，所设计的基于能量反馈控制策略的起摆控制器是有效的，可以实现长摆的自动摆起控制； 长摆经过四次振荡，用时约5．7s便摆到了0°位置附近(即竖直向上位置附近)； 在长摆的起摆过程中，短摆偏离竖直向下位置的最大角度为45°； 当长摆摆到竖直位置附近时，短摆偏离竖直向下位置的角度为10°，处于可控状态。由仿真结果可知，在长摆的起摆过程中，旋转臂偏离初始位置的最大角度为68°； 当长摆摆到竖直位置附近时，旋转臂偏离初始位置的角度为60°，同样处于可控状态。

在以上设计思路中，考虑到短摆的起摆控制是在长摆起摆并稳定控制以后再进行的。因此，在MATLAB/Simulink仿真环境下，搭建了如图3105所示的基于PRDIP可视化模型的仿真模型，以此验证基于能量反馈控制策略的起摆控制器对短摆的控制效果； 图3105中的稳定控制器及其切换条件，将在后文给出。得到的仿真曲线如图3106所示。



图3105短摆的起摆控制仿真模型图




图3106短摆起摆控制仿真结果


由图3106的仿真结果可知，本节所设计的基于能量反馈控制策略的起摆控制器是有效的，可以实现短摆的自动摆起控制； 由于短摆质量较轻，起摆过程所需能量较少，因此短摆仅经过两次振荡，用时约2s便摆到了360°位置附近(即竖直向上位置附近)； 在短摆的起摆过程中，长摆偏离稳定位置的角度在±20°范围内，但是可在稳定控制器的控制作用下，仍然能够保持倒立不倒； 当短摆摆到竖直位置附近时，长摆偏离稳定位置的角度为5°，角度大小符合进入双摆稳定控制的切换条件。进一步，由仿真结果可知，在短摆的起摆过程中，旋转臂偏离初始位置的最大角度为92°，并且当长摆摆到竖直位置附近时，旋转臂偏离初始位置的角度为15°，处于可控状态。

综上，所设计的基于能量反馈控制策略的起摆控制器是有效的，能够实现长摆与短摆的起摆控制，且控制效果符合设计要求。

(2) 基于全状态反馈控制策略的稳摆伺服控制器设计。

系统的稳摆控制即针对PRDIP在竖直不稳定平衡点附近的稳定控制设计，此刻，一般已经对系统的模型进行了线性化，写成状态空间模型形式，即为

x·=Ax+Bu
y=Cx(3102)

式中，状态量为x=［x1,x2,x3,x4,x5,x6］T=［θ,θ·,θ1,θ·1,θ2,θ·2］T
，系统矩阵A和输入矩阵B的表达式参见式(323)，系统输入u=τ
(电机转矩) ，输出选择为三个角度θ、θ1、θ2
。

目前，针对倒立摆的稳摆控制方案有LQR控制(线性二次型最优控制)、滑模变结构控制、模糊控制、H∞控制等方案，但对于平衡点附近线性化模型的稳摆控制，一般采用“全状态反馈”即可得到较好的控制结果； 考虑到系统起摆和稳摆不同状态的控制设计，我们采用LQR控制器实现PRDIP系统的稳摆控制。

LQR控制设计，一般选取如下能量函数作为系统的性能指标

J=12∫∞0(xTQx+Ru2)dt(3103)

式中，Q为与状态误差要求相关的半正定的权重系数矩阵，R为与控制输入能量消耗相关的正定的权重系数(此处由于系统为单输入系统，所以R为标量)。根据最优控制的基本理论，可得到系统的全状态反馈控制器为

u=-Kx=-R-1BTPx(3104)

式中，矩阵P可通过求解以下退化的矩阵Riccati方程得到

PA+ATP-PBR-1BTP+Q=0(3105)

如果存在正定的矩阵P，则表明最优问题式(3103)有解，即可以将求得的正定矩阵P代入式(3104)求得全状态反馈系数K，此时求得的状态反馈系数K即为最优控制。

在实际的控制设计中，可采用MATLAB的lqr函数，结合所获得状态空间表达式(3102)，通过给定最优目标函数(3103)中的权重系数Q和R，即可以求得最佳的状态反馈矩阵K，在仿真分析或实验过程中，可对获得的状态反馈矩阵K中的元素进行微调，以获得最佳的控制效果。

下面，在MATLAB/Simulink仿真平台搭建了如图3107所示的基于PRDIP可视化模型的仿真模型，以此验证所设计的基于LQR的全状态反馈稳摆控制器的有效性。图3107所示的仿真模型中，PRDIP模块为可视化模型封装模块； LongShortStable为基于LQR的全状态反馈控制器。



图3107基于LQR的PRDIP稳摆控制系统仿真图


当给定旋转臂转角的控制目标为θ=0°，长摆的初始角度θ1=0．1rad=5．73°，短摆的初始角度θ2=0．1rad=5．73°时，进行仿真实验，得到的仿真结果如图3108所示。其中，全状态反馈控制器的K矩阵参数取为

K=［-0.11-0.3244.608.69-35.65-3.23］


由图3108可知，在状态反馈控制器的控制作用下，两个摆杆可以稳定在竖直向上位置，且摆角的稳定精度在［-1°,1°］，旋转臂最终也可以稳定在目标位置0°。



图3108稳摆控制仿真结果(初始值为θ=0°，θ1=5．73°，θ2=5．73°)


下面检验所设计的状态反馈控制器能够控制的摆杆的最大初始角度。当给定旋转臂转角的控制目标为θ=0°，长摆的初始角度θ1=0．34rad=20°，短摆的初始角度θ2=0．34rad=20°时，进行仿真实验，仿真结果如图3109所示。



图3109稳摆控制仿真结果(初始值为θ=0°，θ1=20°，θ2=20°)


由图3109可知，当两个摆杆的初始角度增大到20°时，在状态反馈控制器的控制作用下，两个摆杆可以回到并稳定在竖直位置，且摆角的稳定精度在［-1°,1°］； 旋转臂最终也可以稳定在目标位置0°，且角度的稳定精度在［-1°,1°］。但是当一个摆杆的初始角度大于20°时，经过仿真发现，控制器就无法保证PRDIP系统的稳摆控制了。

读者还可以通过设计不同的参数扰动和转矩扰动仿真分析稳摆控制系统的鲁棒性，此处不再单独给出。

(3) “起摆+稳摆”双模态切换全域运动控制器设计。

由前文结合图3102的分析以及相关的起摆和稳摆控制设计可知，串行起摆控制便于将控制器分段； 同时，考虑到长短摆的惯量差距较大，所以采用先起摆长摆，后起摆短摆的方案。设计双模态控制器的切换策略为： 第一控制器为长摆的起摆控制器； 第二控制器为长摆的稳摆控制器； 第三控制器为长摆稳摆控制器与短摆起摆控制器的叠加控制器； 第四控制器为双摆稳摆控制器。同时，为了节约双模态控制的总时长，不妨将第二与第三控制器合并，即稳定长摆的同时，就开始启动短摆的起摆，这也是由于短摆惯量远小于长摆，短摆的起摆转矩需求较小，因此对长摆稳摆造成的干扰较小。

综上，基于对起摆控制器和稳摆控制器的具体设计，以及设定的模态控制控制方案，给出PRDIP系统全域控制器如式(3106)所示，并且给定各阶段控制器的切换角度均为0．2rad。

τout=

-k1(θ·1E1cosθ2+λ1θ·)|θ1|>0.2
-k2sgn(θ·2E2cosθ2+λ2θ·)+K1x|θ1|≤0.2&|θ2|>0.2
K2x|θ1|≤0.2&|θ2|≤0.2

(3106)


为了验证所设计的基于双模态控制器的全域控制方案是否可行，在Simulink平台下搭建了如图3110所示的PRDIP系统全域控制仿真模型。其中控制器表达式(3106)中各参数取值如下： k1=0．29，K1=［1．41 2．16 -34．72 -6．15 0 0］，k2=6．41，K2=［1 2．14 -108．4 -21．09 75．79 6．81］。在仿真模型中，如图3110所示，切换条件的判断通过MATLAB Function函数编写，通过Multiport Switch根据切换条件选择不同的起摆或稳摆控制器。初始值为θ=0°、θ1=180°、θ2=180°的仿真结果如图3111所示。



图3110PRDIP系统全域运动控制仿真模型图




图3111全域运动控制仿真结果曲线


由图3111的仿真结果可知，所设计的基于能量反馈起摆控制/全状态状态反馈稳摆控制的双模态切换全域控制器是有效的，可以实现PRDIP系统从起摆到稳摆的全域控制。初始时，长摆起摆过程大约用时6s； 当长摆角度满足切换条件时，长摆的稳定控制器可将其稳定在竖直位置(θ1=0°)； 当旋转臂回到给定位置附近(θ=0°)且长摆稳定大约2s后，短摆在起摆控制器的作用下，成功实现起摆，起摆过程大约用时2s； 当长摆与短摆的角度同时满足切换条件时，双摆的稳摆控制器可将两个摆杆同时稳定在竖直位置(θ1=0°，θ2=360°)。

进一步可以通过仿真实验检验所设计双模态切换全域控制器的伺服控制能力。当PRDIP系统稳定后，分别在第16s、24s、30s时更改旋转臂目标位置为-45°、-90°、0°，得到的仿真结果如图3112和图3113所示。由图3112的仿真结果可知，PRDIP系统全域控制过程完成后，在全域控制器的控制作用下，旋转臂角度能够跟随给定，在一定时间到达并稳定在目标位置，角度控制精度在［-1°,1°］。由图3113的仿真结果可知，在旋转臂角度跟随3个给定目标位置的过程中，长摆和短摆在全域控制器的控制作用下均能保持直立不倒，长摆的角度波动在［-2°,2°］范围内，短摆的角度波动在［-1°,1°］。因此，通过仿真实验证明了本节所设计的基于能量反馈控制起摆和基于全状态反馈控制稳摆的双模态切换全域控制器对旋转臂位置具有伺服控制能力，能够使旋转臂准确跟随给定位置。



图3112全域位置伺服控制实验旋转臂转角仿真曲线




图3113全域位置伺服控制实验摆杆摆角仿真曲线




拓展阅读
一阶并联旋
转双倒立摆
系统的全域
运动控制
（起摆与稳
摆+三步位
置伺服控制）

4．  结论与思考

本节我们基于MATLAB/Simscape工具箱建立了PRDIP系统的“可视化模型”，应用能量法与LQR(线性二次型最优控制)实现了PRDIP系统“起摆+稳摆+位置伺服”的全域运动控制。有以下几点结论与思考： 

(1) 在3．3节中所建立的PRDIP系统“数学模型”是LQR控制器设计的基础，MBD(基于模型的设计)方法是控制系统设计的基本方法。

(2) 基于MATLAB/Simscape工具箱所建立的PRDIP系统“可视化模型”具有较高的模型精度，可直接用于控制系统(被控对象)的仿真实验； 同时，对于复杂的机电一体化控制系统设计，应用MATLAB/Simscape“可视化模型”实现的仿真实验具有“仿真精度高(建模过程AI自动化)、动态响应直观(运动姿态VR可视化)”的特点，我们应该深入了解与充分应用这一利器。

(3) 本节所述的控制方案是以电机的“转矩”作为系统控制量(控制输入)的，所得到的控制效果(即仿真结果)还是很不错的； 但是，在实物系统装置的控制实现时，在电机驱动器的“控制模式”选择上，我们会面临“转矩控制与转速控制(角速度控制)”两种控制模式中哪一种更合适的问题。在4．4节中，通过深入的模型分析，我们会了解到： 与4．4节所述的把“旋转臂的角加速度”作为PRDIP系统控制量(控制输入)的控制方案相比，后者应该是更好的控制方案； 感兴趣的读者不妨在第4章继续深入研读之。

(4) 本节针对PRDIP系统的全域运动控制问题(起摆+稳摆+位置伺服)仅给出了一种可行的解决方案，实际上还有诸如“串行起摆与并行起摆的时间最优控制问题”“全域运动控制的时间最优控制与能量最优控制问题”，以及串行起摆过程为何要“先长后短”等一些未解问题有待读者回答； 感兴趣的读者可在本节内容的基础上，深入学习与探究之。

读者扫描对应二维码可观看“并行起摆+稳摆+旋转臂位置伺服控制”的实物实验工作视频，进一步思考上述的“时间最优控制”问题。

3.7基于经典频域法的DC/DC变换器控制方案

电力电子变换器是实现电能变换与控制的电力自动化装置，包括DC/DC、AC/DC、DC/AC和AC/AC四种电能变换形式。本节将以Buck变换器为例，讨论基于经典频域法的DC/DC变换器控制系统设计问题。在后续电力电子系统控制问题的讨论中，我们将沿循控制对象从简单的DC/DC变换器到相对复杂的三相电压型PWM整流器，模型特征从单输入单输出系统到多输入多输出系统，控制方法从经典频域法到非线性控制方法的思路展开论述。

1. DC/DC变换器的数学建模

DC/DC变换器可将直流电转换为另一固定电压或可调电压的直流电，有时也称为直流斩波电路，其可应用于开关电源、直流传动系统、电动车能量变换存储等领域。DC/DC变换器包括Buck、Boost、Buck/Boost、Cuk、Sepic和Zeta这6种基本直流斩波电路，以及正激、反激、半桥、全桥和推挽等间接直流变换电路［1419］。本小节将以最常用的Buck型电路为例（如图3114所示），讨论DC/DC变换器的控制系统设计问题。



图3114Buck变换器的拓扑结构


图3114中，Ui表示直流输入电压，uo为直流输出电压，iL为电感电流，R为负载电阻，L和C分别为滤波电感和滤波电容，VD为二极管，V为功率开关器件，s为开关函数，s=1表示功率开关器件开通，s=0表示功率开关器件关断。进而，在连续导电模式（Continuous Conduction Mode，CCM）下，由基尔霍夫电压和电流定律，可建立Buck变换器的数学模型如式(3107)所示。
iL=Cduodt+uoRuo=sUi－LdiLdt(3107)
为便于利用经典控制理论对该系统进行设计，需对式(3107)进行小信号线性化处理。将方程中各变量等效为稳态直流分量和交流小信号扰动值之和的形式，即令uo=Uo+u^o，iL=IL+i^L，s=D+d^。其中，Uo、IL、D为稳态值，u^o、i^L、d^为小信号扰动值。将这些变量代入式(3107)中，并由稳态时各变量的关系，可得微分方程形式描述的系统小信号模型如式(3108)所示，此步骤为“分离扰动”过程。
i^L=Cdu^odt+u^oR
u^o=d^Ui－Ldi^Ldt (3108)
对式(3108)进行拉氏变换，整理后可得Buck变换器的控制输入（占空比）与直流电压输出间的传递函数关系为
Go(s)=Uo(s)D(s)=UiLCs2+LRs+1(3109)
由式(3109)可知，该模型为一双重极点型控制对象，可利用经典控制理论的频域分析方法，设计相应控制器，实现Buck变换器的直流输出电压控制。

2. DC/DC变换器的控制系统设计

（1） Buck变换器被控对象的频率特性。

对于某一Buck型DC/DC变换器，其系统参数为： 直流输入电压Ui=28V，直流输出电压为15V，直流负载电阻R=3Ω，滤波电感L=50μH，滤波电容C=500μF。直流输出电压给定值Uref=1.5V，即电压采样网络H(s)=1/10。对于PWM调制环节Gm(s)=1/Um，其载波信号幅值为Um=1V，开关频率fs=100kHz。在此条件下，基于经典频域法设计该Buck变换器的直流电压控制系统，如图3115所示，其中Gc(s)为待设计的控制器。



图3115Buck变换器控制系统结构



将系统参数代入式(3109)中，利用MATLAB控制系统设计工具箱及相应指令函数，可绘制出被控对象的幅频和相频特性，如图3116所示。



图3116Buck变换器的波特图



系统的开环传递函数为
G(s)=Gc(s)Gm(s)Go(s)H(s)(3110)
当控制器Gc(s)=1时，利用MATLAB命令可绘制系统开环传递函数的幅频特性和相频特性如图3117所示。由该频率特性曲线，可分析Buck变换器的稳定性、稳态性能和动态性能。



图3117Gc(s)=1时开环传递函数的波特图

首先，由图3117可知，系统的剪切频率为ωc=12.3krad/s，相角裕度为φm=4.2°。相角裕度较低（接近于零），使得系统虽然理论上是稳定的，但当Buck变换器受到一定的参数摄动或外部扰动时，系统将会变得不稳定； 其次，系统的直流增益Gu0=HUi/Um=2.8，据此可计算出稳态误差为1/1+Gu0=26.3%，如此大的稳态误差是不能满足实际应用要求的； 最后，系统的剪切频率为ωc=12.3krad/s，对应fc=1.96kHz，剪切频率过小，使得系统的动态响应速度很慢。综上所述，未施加有效控制时，Buck变换器不能满足实际应用中“稳、准、快”的需求，需要设计合理的控制器，提高系统的综合性能。



（2） 电力电子变换器频域法设计的原理与步骤。

在利用经典频域法设计电力电子变换器控制系统时，主要用波特图来表示被控对象、控制器以及开环传递函数的频率特性。系统开环传递函数的波特图可以反映出闭环系统的稳定性、稳定裕度、稳态和动态性能。理想的开环传递函数频率特性在低频段、中频段和高频段应该分别满足如下要求［1921］： 

① 低频段： 开环传递函数频率特性的低频段反映了系统包含积分环节的个数和直流增益的大小，主要影响系统的稳态性能。对于DC/DC变换器等电力电子系统，理想的低频特性是直流增益无限大，低频段以-20dB/dec的斜率下降。

② 中频段： 开环传递函数频率特性的中频段需以-20dB/dec斜率下降并穿越0dB线，剪切频率（或称穿越频率）与系统的稳定性、调节时间和超调量等动态性能密切相关。

③ 高频段： 高频段与系统的稳态和动态性能关系不大，但其反映了系统对高频干扰信号的抑制能力。高频段幅频特性衰减越快，系统的抗干扰能力就越强，一般要求以-40dB/dec斜率下降。

电力电子变换器频域法设计的一般步骤可概括为： 把系统的性能指标和技术要求转换为开环传递函数的波特图，即期望特性； 根据被控对象的波特图和开环传递函数的波特图确定控制器的波特图； 基于控制器的波特图，选择合适的控制器并完成参数设计。

下面，基于电力电子变换器频域法设计的原理与步骤，给出Buck型DC/DC变换器控制系统的设计过程。

（3） Buck变换器的PD控制器设计。

针对前述Buck变换器在控制性能方面的不足，这里设计合理的控制器进行改进。提高系统相角裕度的一个有效办法是采取PD控制器（超前校正装置）。此时，在小于系统剪切频率处，给控制器增加一个零点，使开环传递函数产生足够的超前相移，可保证系统获得较大的相角裕度。另一方面，在大于剪切频率处，给控制器增加一个极点，提高开环传递函数高频段的下降斜率，可更好地抑制高频噪声，同时这种控制器也便于工程实现。

PD控制器的传递函数如下式(3111)所示
Gc(s)=Gc01+s/ωz1+s/ωp(3111)
式中，ωz<ωc<ωp。

需要说明的是，式(3111)所示的PD控制器在文献［1921］中也称为PD补偿网络或超前校正装置； 本节所说的“控制器”与“补偿网络/校正装置”具有同一含义。实际上，由于理想微分环节不可物理实现，式(3111)所示的PD控制器与传统意义上的理想比例微分控制器稍有差别； 再者，由于ωz<ωp，系统仅在中低频段才具有理想的比例微分特性。

为了提高剪切频率，设加入PD控制器后，开环传递函数新的剪切频率f′c为开关频率fs的1/20，即f′c=fs/20=5kHz。设补偿后新的相角裕度φ′m=52°，则PD控制器的零点、极点以及直流增益计算公式为［1921］

ωz=ω′c1－sinφ′m1+sinφ′m=2π×5×1－sin52°1+sin52°=10.8k rad/s(3112)
ωp=ω′c1+sinφ′m1－sinφ′m=2π×5×1+sin521－sin52=91.2k rad/s(3113)
Gc0=ωzωpLC(ω′c)2/Gu0
=10.891.2×50×10－6×500×10－6×(2π×5×103)2/2.8=3(3114)

此时系统的开环传递函数为
G(s)=Gc(s)1UMGo(s)H(s)
=Gc0UiHUM1+sωz1+sωpLCs2+LRs+1(3115)

利用MATLAB命令可分别绘制PD控制器和系统开环传递函数的频率特性曲线如图3118和图3119所示。



图3118PD控制器的波特图




图3119采用PD控制器后系统开环传递函数的波特图



由图3118和图3119可知，采用PD控制器后，实际系统的相角裕度为φ′m=53.2°，剪切频率为ω′c=32.2k rad/s，即f′c=5.12kHz。同时对于中频段的10k rad/s 至100k rad/s较大范围，相角裕度均维持在40°以上，可保证系统受到较大参数摄动和外部扰动时，仍然具有较好的稳定性和动态性能。

然而，系统幅频特性曲线在低频段较为平直，会存在较大的稳态误差，有必要进一步对PD控制器进行改进。

（4） Buck变换器的PID控制器设计。

为了解决PD控制器存在的系统稳态误差较大的问题，可在其传递函数基础上，通过加入倒置零点，改善开环传递函数的低频特性，构成如式(3116)所示的PID控制器。
Gc(s)=Gc11+s/ωz1+ωm/s1+s/ωp(3116)
与前述PD控制器相比，PID控制器仅在低频段有所改变，而在中高频段特性不变，因此所设计的倒置零点的频率应远离系统的剪切频率，这样系统的相角裕度和剪切频率可基本维持原值，不受其影响； 控制器的零点ωz、极点ωp以及直流增益Gc1也可沿用原有值。倒置零点的频率一般取为
ωm=ω′c/10=3.22k rad/s(3117)
此时系统的开环传递函数为
G(s)=Gc(s)1UMGo(s)H(s)
=Gc0UiHUM1+sωz1+ωms1+sωpLCs2+LRs+1(3118)
利用MATLAB命令可分别绘制PID控制器和系统开环传递函数的频率特性曲线如图3120和图3121所示。



图3120PID控制器的波特图




图3121采用PID控制器后系统开环传递函数的波特图


由图3120和图3121可知，采用PID控制器后，系统的中高频特性基本保持不变，系统的相角裕度为φ′m=47.5°，剪切频率为ω′c=32.3k rad/s，即f′c=5.14kHz。 在低频段，系统开环传递函数可近似为式(3119)所示的积分环节，其幅频特性曲线以-20dB/dec的斜率下降，保证了系统的稳态性能。至此完成了Buck变换器控制系统的设计。
G(s)≈Gc0UiHωmUM1s(3119)

3. 仿真实验

根据Buck变换器系统参数及前述设计的控制器传递函数，可在MATLAB/Simulink环境下，建立Buck变换器控制系统仿真模型，如图3122所示。该仿真模型包括利用SimPowerSystems仿真工具箱搭建的Buck变换器主电路部分，以及由电压给定环节、控制器、PWM调制环节和检测环节组成的控制电路部分。利用该仿真模型，可对所设计的Buck变换器控制系统进行仿真实验验证。实际上，本节第二部分所绘制的系统波特图即是利用MATLAB指令来获取的，但其侧重于系统的频域分析，本部分将给出系统时域下的仿真结果，其与频域分析存在着一定的对应关系。



图3122基于PID控制器的Buck变换器控制系统仿真模型


仿真模型中，由于反馈电压采样网络H(s)=1/10，因此Buck变换器的输出侧若期望获得15V的直流电压，则控制系统给定值Uref=1.5V，图3123为采用PID控制器和PD控制器的直流输出电压响应曲线。




图3123Buck变换器的直流电压输出曲线


由仿真结果可以看出，采用PID控制器后，Buck变换器的直流输出电压在t=0.3ms即达到期望值，具有较快的响应速度，同时无稳态静差。然而，采用PD控制器，虽然Buck变换器可实现稳定控制，并具有较好的快速性，但稳态时直流输出电压仅为13.34V，存在着1.66V的稳态静差。这一实验结果与此前的频域分析是相符的。

需补充说明的是，仿真模型中滤波电容的电压初始值设为9V，以防止初始上电状态下，由于电容电压不能突变而导致的直流输出电压较大过冲，这与实际装置中一般需配置的电容预充电环节是相符的。

为了更好地验证所设计的基于PID控制器的Buck变换器控制系统性能，这里考虑系统受到直流负载扰动和直流输入电压扰动两种情况，相关仿真实验结果如图3124所示。



图3124Buck变换器控制系统的抗扰性能



由仿真结果可见，t=0.6ms时，直流负载从R=3Ω变化为R=1.5Ω，此时仿真实验结果如图3124(a)所示。t=0.6ms时，直流负载保持R=3Ω不变，直流输入电压从Ui=28V变化为Ui=21V，此时仿真实验结果如图3124(b)所示。可见，对于直流负载和直流输入电压大范围的扰动，Buck变换器都能较为快速地回到稳定状态，系统具有较好的抗扰能力。

4. 小结

本节以Buck变换器为例，给出了基于经典频域法的DC/DC变换器控制系统设计过程，综合起来有以下几点结论： 

（1） 本节给出了基于经典频域分析理论的电力电子变换器控制系统设计过程，总结了理想开环传递函数的频率特性（低频段、中频段和高频段）与控制系统稳定性、相对稳定性（相角裕度）、稳态性能（稳态误差）和动态性能（快速性等）间的关系，为电力电子变换器控制系统的频域分析与综合提供了理论指导。

（2） 以波特图为工具，对Buck型DC/DC变换器这类单输入单输出系统，利用经典频域法设计其控制器，并进行了仿真实验验证。仿真结果表明，一方面，系统的频域分析与时域响应可互为验证； 另一方面，所设计的Buck变换器控制系统具有较好的稳态、动态和抗扰性能。

（3） 本小节设计了Buck变换器的电压单环控制系统，单环系统的优点在于结构简单、设计方便，但也存在扰动下系统响应速度慢等问题（例如图3124(b)所示直流输入电压扰动后的输出电压曲线）。实际上，为了追求更好的控制性能，在直流输出电压环基础上，可引入电感电流进行内环控制（包括峰值电流控制、平均电流控制和滞环电流控制等），构成双环控制系统。此时，控制系统的设计要更复杂些。

（4） 在进行DC/DC变换器控制系统实物设计时，还需考虑控制器的硬件实现问题，因其不属于本书的讨论范围，此方面未过多展开论述，这里只简要给出其两种实现途径： 一种方法是“模拟实现”，即采用有源校正装置，利用运算放大器和电阻、电容来实现； 另一种方法是“数字实现”，近年来随着高性能处理器的飞速发展，这种方法已成为技术主流。

需要说明的是，本节所阐述的Buck变换器控制系统设计并未考虑滤波电容的串联等效电阻，以及不连续导电模式（Discontinuous Conduction Mode，DCM）工况，对此感兴趣的读者可基于本节内容进一步自行设计完成。

3.8三相电压型PWM整流器的高功率因数控制方案

整流器作为电力电子设备的前端电路，应用极其广泛。传统的整流装置采用二极管不控整流或晶闸管相控整流方式，具有网侧电流谐波大、功率因数低等缺点，给电网带来了严重的谐波和无功功率“污染”问题。三相电压型PWM整流器（简称“PWM整流器”）采用全控型电力电子器件，利用PWM斩波控制方式，能够实现交流侧高功率因数运行，具有交流电流畸变小、输出直流电压可调以及能量可双向流动等优点，越来越受到人们的广泛关注［2226］。

PWM整流器稳态和动态控制性能的优劣依赖于控制系统的合理设计，PWM整流器控制系统设计包括“控制系统结构设计”和“控制器/策略设计”两部分； 在控制系统结构设计方面，本节基于欠驱动系统理论，给出了PWM整流器“有功/无功电流内环、直流电压外环”的双闭环控制系统结构，并阐述了该结构设计的理论成因； 在控制器/策略设计方面，由于PWM整流器是非线性系统，采用线性系统理论设计的控制策略存在着“控制器参数整定困难、大范围扰动不稳定”等问题，故本节将采用非线性控制方法——“滑模变结构控制方案”，以实现PWM整流器的高功率因数控制。

1. 系统建模与模型验证

在第2章2.5.6节“PWM整流器控制问题”中，我们给出了三相静止坐标系下PWM整流器的数学模型，这里我们用“必要条件法”对所建模型进行验证。所谓必要条件法，就是所进行的模型验证实验的结果是依据经验可以判定的，其正确性的结果是正确的模型所应具备的必要性质。

在MATLAB/Simulink环境下，按照PWM整流器数学模型，搭建其仿真模型如图3125所示，这里函数Fcn，Fcn1，Fcn2，Fcn3和Fcn4的表达式分别为

Fcn： （u［4］-u［11］*u［1］-u［10］*u［7］+u［10］*(u［7］+u［8］+u［9］)/3）/u［12］

Fcn1： （u［5］-u［11］*u［2］-u［10］*u［8］+u［10］*(u［7］+u［8］+u［9］)/3）/u［12］

Fcn2： （u［6］-u［11］*u［3］-u［10］*u［9］+u［10］*(u［7］ +u［8］+u［9］)/3）/u［12］

Fcn3： （u［1］*u［7］+u［2］*u［8］+u［3］*u［9］-u［10］ /u［14］）/u［13］

Fcn4： u［1］*u［7］+u［2］*u［8］+u［3］*u［9］



图3125模型验证仿真模型



系统模型方程中的常量R,L,C,RL是可变参数，其可在模型外部进行灵活设置，这里取R=0.3Ω，L=20mH，C=990μF，RL=200Ω。由idc=iasa+ibsb+icsc，可得不同开关模式时的idc值（如表37所示）。


表37PWM整流器不同开关模式时的idc取值



scsbsa
001
010
011
100
101
110
111
000
idc
ia
ib
－ic
ic
－ib
－ia
0
0


为了检验PWM整流器数学模型是否与实际系统相符，下面设计了3个仿真实验，如表38所示。


表38模型验证仿真实验



序号
实验条件
判定依据
实验一
sa=sb=sc=0
idc=0
实验二
sa=1,sb=sc=0
idc=ia
实验三
sa=0,sb=sc=1
idc=－ia


模型验证的仿真实验结果如图3126所示。由实验结果可以看出，该模型行为与理论分析完全符合，因而可在一定程度上证明所建模型的正确性。



图3126模型验证的仿真实验结果


2. PWM整流器的双闭环控制系统构成

（1） PWM整流器的欠驱动特性。

在第2章2.5.6节“PWM整流器控制问题”中，给出了dq同步旋转坐标系下PWM整流器的数学模型，如下式所示： 
Ldiddt=－Rid+ωLiq－sduo+ed
Ldiqdt=－ωLid－Riq－squo+eq
Cduodt=sdid+sqiq－uoRL(3120)
由该模型可知，描述其位形空间内运动的独立变量为有功电流id、无功电流iq和直流电压uo，因此系统的自由度为3。另一方面，PWM整流器的控制输入有sd和sq 2个。由于其控制输入个数小于系统自由度，因此PWM整流器是一个“欠驱动电力电子系统”，这样可在欠驱动系统的理论体系下，对PWM整流器进行重新审视与讨论。

“欠驱动”这一概念最早用于描述机械系统的运动控制问题，所谓“欠驱动系统”指的是系统控制输入数目小于自由度的系统； 例如，对于一阶直线倒立摆系统，其控制输入只有1个，即小车水平方向控制力，但系统的自由度为2（小车位移和摆杆摆角）； 实际上，龙门吊车、自平衡两轮电动车、水面舰船等系统都具有同类性质； 这种控制输入的缺失，给“欠驱动系统”的有效控制带来困难。

（2） 双闭环控制系统结构设计与分析。

对于欠驱动机械系统，驱动变量与欠驱动变量的选择可根据系统运动方式直接判断得出。例如，对于一阶倒立摆系统，小车位移为驱动变量，摆杆摆角为欠驱动变量。此处，“驱动变量”指的是在控制输入激励下，可以直接实现稳定控制的状态变量； 所谓“欠驱动变量”指的是控制输入并不直接激励，而仅由系统内部动态(零动态)影响的状态变量。

然而，欠驱动电力电子系统的驱动变量与欠驱动变量的选择则需要对系统稳定性进行分析后方能确定。亦即，对于PWM整流器，由于其仅有2个控制输入，因此状态变量包括2个驱动变量(与sd、sq这2个控制输入对应)和1个欠驱动变量(无控制输入直接激励)。

选取id、iq、uo中哪2个状态变量作为驱动变量需进行理论分析，这里有3种选择方案： 

① 选择uo、iq为驱动变量，id为欠驱动变量。

由于采用了等功率坐标变换，因此PWM整流器的交直流侧瞬时功率平衡关系式为
edid=Liddiddt+iqdiqdt+R(i2d+i2q)+Cuoduodt+u2oRL(3121)
式(3121)实际上反映了PWM整流器能量变换的本质特征。由该式可分别获得id、iq、uo的内部动态(零动态)方程，进一步可判断出欠驱动零动态子系统，乃至整个PWM整流器控制系统的稳定性［2728］。

PWM整流器的控制目标，一是实现uo收敛于给定值u*o(即直流电压稳定并可调)，二是实现iq收敛于给定值i*q=0(即实现网侧单位功率因数控制)。显然，选择uo、iq作为驱动变量，可直接实现这2个目标，从而欠驱动变量id的稳定性就直接决定PWM整流器控制系统的稳定性。

在式(3121)中，令iq=0，uo=u*o，可得以id为欠驱动变量的系统零动态方程为
diddt=－RidL－(u*o)2LidRL+edL(3122)


图3127欠驱动变量id的相平面图

此时，可得如图3127所示的id零动态相平面图。

此时，系统有Id+和Id- 两个平衡点，其中Id+是稳定的平衡点，Id-是不稳定的平衡点。Id+的表达式为
Id+=12edR+e2dR2－4(u*o)2RRL(3123)
由式(3123)可见，Id+远远超出了PWM整流器的运行范围，物理上不可实现，因此欠驱动变量id只存在唯一的不稳定平衡点Id-。由于id零动态不稳定，故选取id作为欠驱动变量是不可行的。

② 选择uo、id为驱动变量，iq为欠驱动变量。

在该方案中，施加有效控制后，uo和id将分别达到给定值u*o和i*d，从而欠驱动变量iq的稳定性就直接决定PWM整流器控制系统的稳定性。

在式(3121)中，令uo=u*o，id=i*d，则可得以iq为欠驱动变量的系统零动态方程： 

edi*d=Liqdiqdt+R((i*d)2+i2q)+(u*o)2RL(3124)
由稳态时功率平衡关系
edi*d=R(i*d)2+(u*o)2RL(3125)
可得以iq为欠驱动变量的系统零动态方程为
diqdt=－RLiq(3126)


图3128欠驱动变量iq的相平面图

此时，iq的零动态相平面图如图3128所示。图中，Iq=0是无功电流iq的稳定平衡点，因此理论上可只对id、uo进行直接控制，间接使得系统内部动态iq收敛到零。然而实际系统往往存在参数不确定性、未知扰动和未建模动态，使得iq自身归于零的收敛过程存在稳态误差，难以获得较高的功率因数。因此对于我们更加关注的无功电流iq的控制问题，采用直接控制更为合适。


另一方面，即使想将iq置于外环，通过i*d或u*o进行间接控制，但由式(3127)所示的id和iq间传递函数关系式可知，由于i*q=0，即式(3127)分母为零，故无法通过调节i*d来确保iq得到有效控制。因此，选取iq作为欠驱动变量也是不可行的。
Iq(s)=ed－2Ri*d－sLi*d(2R+sL)i*qId(s)(3127)
③ 选择id、iq为驱动变量，uo为欠驱动变量。

在该方案中，sd、sq施加有效控制后，id和iq将分别收敛于给定值i*d和i*q=0。进而欠驱动变量uo的稳定性就直接决定PWM整流器系统的稳定性。

在式(3121)中，令id=i*d，iq=0，则可得以uo为欠驱动变量的系统零动态方程为
duodt=1Cuoedi*d－R(i*d)2－u2oRL=1RLCuo((u*o)2－u2o)(3128)


图3129欠驱动变量uo的相平面图

其相平面图如图3129所示。由该图可看出，Uo是直流侧电压的稳定平衡点(由于-Uo<0，直流电压不可能为负，故这里不必考虑-Uo)。因此，理论上可只对id、iq进行直接控制，间接使得uo收敛到稳态值u*o。同样，为了解决实际系统uo存在稳态误差及收敛速度不可控的问题，需要引入直流电压uo的外环控制器，通过对该控制器输出i*d的调节，间接实现uo迅速收敛到u*o。

对PWM整流器模型进行小信号处理，可得id和uo间传递函数关系式如式(3129)所示。由此可见，通过设计电压外环控制器，并将其输出作为id的给定值i*d，进而通过内环i*d的调节，间接实现uo的有效控制是可行的。
Uo(s)=RL(ed－2Ri*d－sLi*d)(2+sRLC)u*oId(s)(3129)
由以上分析可知，所确定的PWM整流器双闭环控制系统结构为： 选择id、iq作为驱动变量，利用sd、sq进行内环直接控制； 选择uo作为欠驱动变量，设计外环控制器，利用i*d对其进行间接控制(如图3130所示)。



图3130PWM整流器的双闭环控制系统结构


3. PWM整流器的滑模变结构控制器设计

在已确定的PWM整流器控制系统结构基础上，这里进一步讨论PWM整流器的滑模变结构控制器设计。由于在“龙门吊车重物防摆的滑模变结构控制方案”一节中，已对滑模变结构控制的相关概念和理论进行了介绍，因此本节不再对此赘述，而是重点探讨滑模变结构控制在电力电子系统中的应用与实现问题。

（1） 电流内环控制器的设计。

在电流内环控制器设计时，暂不考虑式(3120)中的直流电压方程，即此时PWM整流器的电流环数学模型为
diddt=1Led－1Lsduo+ωiq－RLid
diqdt=1Leq－1Lsquo－ωid－RLiq(3130)
在此数学模型基础上，我们将设计滑模变结构控制器，使id和iq都能快速达到给定值。

① 定义两个滑模面： 
s1=α1(id－i*d)
s2=α2(iq－i*q)(3131)
其中，α1,α2均是正实数，i*d,i*q分别是id,iq的给定值。

② 采用等速趋近率： 
dsdt=－εsgn(s)(3132)

其中，ε为正实数，则可保证s·s·<0，满足广义滑模条件： 
s·1=－W1sgn(s1)
s·2=－W2sgn(s2)W1>0,W2>0(3133)
③ 求出控制律： 

因为
s·1=α1i·d=α11Led－1Luosd+ωiq－RLid
s·2=α2i·q=α21Leq－1Luosq－ωid－RLiq(3134)
所以由式（3133）和式（3134），可得电流内环控制器的控制律如下： 
sd=Luo1Led－RLid+ωiq+W1α1sgn(s1)
sq=Luo1Leq－RLiq－ωid+W2α2sgn(s2)(3135)
（2） 电流内环控制器的改进。

由式（3135）可以看出，电流内环控制器的控制律存在不连续符号函数sgn(s)，与龙门吊车系统的滑模变结构控制一样，sgn(s)的存在将带来滑模变结构控制所特有的“抖振”问题。

这里，我们采用“边界层法”，即用连续的饱和函数sat(s)来代替不连续函数sgn(s)。饱和函数sat(s)可写为如下形式： 
sat(s)=+1s>Δ
kss≤Δ
－1s<－Δ(3136)
式中，Δ被称为“边界层”，且Δ·k=1。sgn(s)和sat(s)可用图3131所示曲线表示。 



图3131sgn函数及其饱和函数连续化曲线


因此，改进后的电流内环控制器如下所示： 
sd=Luo1Led－RLid+ωiq+W1α1sat(s1)(3137)
sq=Luo1Leq－RLiq－ωid+W2α2sat(s2)(3138)
（3） 电压外环控制器的设计。

PWM整流器控制系统的电压外环通常采用PI控制算法，使得输出直流电压无稳态误差。由于控制系统采用双闭环结构，因此根据电流内环的闭环传递函数，即可设计电压外环PI控制器。可将电流内环和PWM调制环节等效为一个一阶惯性环节［29］： 
Wci(s)=0.753Tss+1(3139)
其中，Ts为电流内环采样周期，也就是PWM开关周期。电压外环控制系统结构如图3132所示。



图3132电压外环控制系统结构


这里，电压外环控制器为
Wpi(s)=Kpiτis+1τis(3140)
则系统开环传递函数为
G(s)=0.75Kpi(τis+1)Cτis2(3Tss+1)(3141)
由此，得电压环中频宽h为
h=τi3Ts(3142)
按照“典型Ⅱ型系统”控制器参数整定方法［1］，可得
0.75KpiCτi=h+12h2(3Ts)2(3143)
选择h=5，可以得到电压外环的PI调节器系统参数如下： 
τi=5×3Ts(3144)
Kpi=4C3Ts(3145)

在实际工程应用中，需根据上述工程化设计方法，确定电压外环PI调节器的基本参数； 在此基础上通过适当的参数调整，就可以找到一组系统稳态/动态性能指标兼顾的PI参数，这就减少了实际工程中PI调节器参数整定的盲目性。

综上，我们得到如图3133所示的PWM整流器滑模变结构控制系统结构，下面对其进行仿真实验。



图3133PWM整流器滑模变结构控制系统结构图


4. 仿真实验

（1） PWM整流器滑模变结构控制系统的仿真。

根据SimPowerSystems仿真工具箱，可在Simulink环境下建立PWM整流器的仿真模型如图3134所示。这一模型将作为被控对象，在控制系统仿真时被使用。



图3134PWM整流器的仿真模型


根据前面设计的控制器表达式，在MATLAB/Simulink环境下，建立PWM整流器滑模变结构控制系统的仿真模型，如图3135所示。



图3135PWM整流器滑模变结构控制系统仿真模型



该仿真模型中，PI controller为电压外环PI控制器，为保证系统稳定运行，对其输出进行了饱和限幅处理。dq_to_abc Transformation为控制信号从同步旋转坐标系到三相静止坐标系的“等功率”坐标变换矩阵。abc_to_dq Transformation1和abc_to_dq Transformation2分别为三相电流和三相电压从三相静止坐标系到同步旋转坐标系的“等功率”坐标变换矩阵。SPWM为PWM调制模块，PWM Converter为被控对象（PWM整流器）的仿真模型，是将图3134所示模型进一步封装得到的。Fcn1和Fcn2为式（3137）和式（3138）所示的电流内环控制器，其值分别为： 

Fcn1： u［12］*(u［1］/u［12］-u［13］/u［12］*u［3］+u［14］*u［4］+u［10］*u［9］/u［8］)/u［5］

Fcn2： u［12］*(u［2］/u［12］-u［13］*u［4］/u［12］-u［14］*u［3］+u［11］*u［7］/u［6］)/u［5］

PWM整流器参数为R=0.3Ω，L=0.02H，C=990μF，ea=78cosωt，eb=78cos(ωt－120°)，ec=78cos(ωt+120°)，负载电阻RL=100Ω。对于滑模变结构控制器，其参数α1=0.765,α2=0.5,W1=15000，W2=4500，k1=k2=14.5。对于PI控制器，给定电压u*o=250V，控制器参数Kpi=0.188，τi=0.0687，饱和限幅值Vsat=±10V。需要注意的是，这里的PI参数是在式（3144）、式（3145）的理论计算基础上，进一步通过仿真实验试凑得到的最佳参数。控制系统仿真实验结果如图3136所示。



图3136PWM整流器滑模变结构控制系统的仿真曲线



从上述仿真结果可以看出： (1)系统具有较快的动态响应速度，有功、无功电流和直流侧输出电压于t=0.08s时达到稳态期望值； (2)从直流侧输出电压波形可以看出，稳态时输出直流电压无静差； (3)从功率因数波形可以看出，在t=0.01s之前，由于此时存在电流波形畸变，导致功率因数较低（在0.966左右），达到稳态之后，功率因数为0.9994，实现了单位功率因数控制。

为更好地验证PWM整流器控制系统性能，这里进一步进行负载扰动实验和给定电压突变实验。在t=0.15s时，负载由RL=100Ω突变为RL=87.5Ω（即再并联一个700Ω的电阻）； 在t=0.45s时，直流给定电压由u*o=250V突变为u*o=280V，仿真实验结果如图3137所示。从仿真实验结果可以看出，所设计的PWM整流器滑模变结构控制系统对于负载变化具有很好的抗扰性能，同时对于给定直流输出电压的变化表现出了很快的动态响应速度。



图3137负载扰动及给定电压突变实验



（2） 控制方案比较分析。

对于上述实验结果，我们与传统的“双闭环PID控制方案”［3032］相比较，可以得出以下几点结论： 

① 滑模变结构控制系统的稳态/动态性能指标要优于传统PID控制。

由于滑模变结构控制的滑动模态可以自行设计，并且对于系统参数及扰动具有不变性，这就使得滑模变结构控制具有快速动态响应、物理实现简单、对参数变化及扰动不灵敏、鲁棒性强等优点。对于PWM整流器，这种强鲁棒性表现为对于负载扰动和电网电压扰动具有很强的抗扰性能。

② 与传统PID控制相比，滑模变结构控制器的参数整定更易于实现。

滑模变结构控制属于非线性控制策略，可实现PWM整流器系统大范围工作稳定，因此其参数整定较为容易。而传统PID控制一般基于系统小信号线性化模型进行设计，PID控制器参数仅仅局限于稳态工作点。

当然，滑模变结构控制也存在一些缺点，最突出的问题就是控制力的“抖振”现象。在滑模变结构控制的具体工程应用中，我们可以采用如前所述的一些“消抖”方法，使抖振幅度限定到工程允许范围内。然而，“抖振”现象只能在一定程度上抑制，而无法从根本上消除； 消除了抖振，也就消除了滑模变结构控制所特有的强鲁棒性。

5. 小结

本节针对PWM整流器的高功率因数控制系统进行了设计与仿真实验，综合起来有以下几点结论： 

（1） 以三相电压型PWM整流器为例，探讨了电力电子系统的建模、模型验证、控制器设计及控制系统仿真等问题； 应用“SimPowerSystems仿真工具箱”可以较为方便地完成电力电子变换器控制系统仿真，加快电力电子系统研究/开发进度。

（2） 讨论了PWM整流器的“欠驱动特性”，并借助欠驱动系统的概念与分析方法，结合零动态稳定性分析基本原理，对PWM整流器双闭环控制系统结构的设计成因进行了阐述与诠释。

（3） 设计了基于滑模变结构控制的PWM整流器控制系统，并进行了仿真实验研究。仿真结果表明，相比于传统PID控制方法，本节所设计的滑模变结构控制系统具有较快的动态响应速度、较强的抗负载扰动能力和给定电压突变跟随性能，可实现高功率因数运行。同时，该滑模变结构控制系统易于控制器参数整定，便于工程实现。

需要说明的是，本节所阐述的PWM整流器控制系统设计是在假定电网电压平衡的条件下完成的。然而，实际工况中三相电网电压在幅值和相位上往往是不平衡的。因而，研究“电网不平衡条件下”PWM整流器的高功率因数控制问题具有重要的实际意义，同时也具有很大的挑战性。感兴趣的读者可参阅文献［3337］对此问题展开研究。

3.9问题与探究——灵长类仿生机器人运动控制研究
1. 问题提出

人类文明的历史也就是人类认识、利用和改造包括自身在内的自然的历史。在这个过程中,由于自身能力的局限性，人类首先间接地延长了个人的肢体，制造了从石器到青铜器、铁器等各种末端执行工具； 其次间接地延长了人类的感觉器官，制造出了具有视觉、听觉、触觉、味觉、嗅觉等更加强大的各种感知工具； 
又进一步间接地延长了人类的大脑，制造出了以电子计算机为代表的各种信息处理和计算工具。人类能力的这种间接延长在20世纪最有说服力的成就
也就是当代最高意义上的自动化的产物——机器人。而人类能力的间接延长的最高产物或者说


图3138长臂猿Ⅰ

21世纪的机器人的发展阶段就是“仿生机器人”。


仿生机器人是仿生学的各种先进技术与机器人领域的各种应用目的的最佳结合。从机器人的角度来看，仿生机器人是机器人发展的最高阶段； 从仿生学的角度来看，仿生机器人是仿生学技术的完美综合与全面应用。仿生机器人既是机器人研究的最初目的,也是机器人发展的最终目标之一。

仿生机器人技术是当今机器人与运动控制的前沿课题，而灵长类(如图3138所示)仿生机器人的运动控制(如图3139所示)又是这一课题中的难点。



图3139长臂猿的运动特点


目前，在灵长类仿生机器人的研究领域，日本名古屋大学的福田敏男教授及其课题组处于领先位置。图3140~图3142所示为其研制的双摆式仿生长臂猿仿生机器人。图3140所示机器人具有两个自由度、一个驱动关节，并在每杆前段带有由电机驱动的手爪和夹子。

 通过对长臂猿生理特性的研究，可以发现，其运动方式主要是以双臂悬垂于树枝上通过自身的荡跃而前进。基于此可以将其简化为一两杆相连的机器人模型(如图3143所示)。

通过简化后的机构模型可以看出，该仿生机器人系统具有两个自由度但只有一个驱动关节。因此，这也是一种典型的“欠驱动机械系统”。




图3140长臂猿Ⅱ




图3141有毛皮的长臂猿Ⅲ




  图3142脱去毛皮的长臂猿Ⅳ





2.　系统建模


基于图3143所示的机构原理，可有图3144所示的系统动力学原理，利用分析力学中的拉格朗日方程，可推导出系统的运动方程。




图3143仿生长臂猿的简化机构模型




图3144系统动力学建模原理图




根据拉格朗日方程,可以求得机器人系统的动力学方程为

M(q)q¨+C(q¨,q)q¨+g(q)=Sτλ(3146)


式中，q,q¨,q¨∈Rn分别是关节点的位移、速度和加速度； M(q)是n ×n 惯性矩阵； C(q·,q)
表示哥氏力和离心力项； g(q)表示重力项； τλ∈Rn表示作用于系统关节的力矩,即系统控制输入。S为输入矩阵,若S可逆，则系统是全驱动的,反之是欠驱动的。在欠驱动情况下，矢量q可分割为q=(q1,q2)，其中q1∈Rn-m，q2∈Rm分别表示被驱动状态和驱动状态。对
M(q),C(q¨,q),g(q)进行相应的分割,可得如下运动方程： 


M11q¨1+M12q¨2+h1+1=0
M21q¨1+M22q¨2+h2+2=τλ(3147)

在忽略系统的关节摩擦力和随机干扰时,其动力学方程可表示为

m11(q)q¨1+m12(q)q¨2+h1(q,q·)+g1(q)=0
m21(q)q¨1+m22(q)q¨2+h2(q,q·)+g2(q)=τ
(3148)

其中

m11(q)=m1L2g1+m2L21+m2L2g2-2m2Lg2L1 cosq2+
I1+I2(a)
m22(q)=m2L2g2+I1(b)
m12(q)=m21(q)=(m2L2g2-m2Lg2L1cosq2)+I2(c)
h1(q，q·)=m2Lg2L1q·2（2q·1+q·2) sinq2(d)
h2(q，q·)=-m2(q·1
)2Lg2L1 sinq2 (e)
g1(q)=(m1gLg1-m2gL1)cosq1-m2gLg2cos(q1+q2)(f)
g2(q)=-m2gLg2cos(q1+q2)(g)


3. 问题探究

(1) 建模问题。

式(3148)给出的是系统模型的微分方程形式，为了便于人们运用已有的控制理论与方法来进行控制系统的设计，往往还需要通过简化与近似，建立起系统的状态方程模型(如式(3149))与传递函数模型(如图3145所示)。

x·=Ax+Bu
Y=Cx+Du
(3149)



图3145系统的传递函数模型



式(3149)中A、B、C、D矩阵为何？图3145中方框内的G1(s)与G2(s)又为何？请读者来回答一下。


(2) 控制问题。

系统的控制问题可由图3146所示的机理来描述，为分析方便将图3144中的q1,q2重新定义为θ1=360°-q1,θ2=q2
。这里有如下两个问题： 



图3146控制问题的机理



Acrobot问题当系统通过适当的控制，使系统状态满足θ1∈(0-δ,0+δ)θ2∈(0-δ,0+δ),δ→0时,该系统属于一类“Acrobot问题”的控制问题。


运动控制问题当系统通过适当的控制，使系统状态满足
θ1∈(A-δ,A+δ)

θ2∈(B-δ,B+δ),δ→0且
θ·1→0
θ·2→0时,我们认为该仿生机器人的末端可有效抵达“预置目标”，即，该问题属于一类“可达目标”的运动控制问题。

针对上述问题，读者不妨应用(或自学)“人工智能控制”与“非线性控制”等理论与方法进行探究。

本章小结

本章在前述内容的基础上，应用MATLAB/Simulink软件及其工具箱，在如下几方面开展了数字仿真技术的应用研究： 

(1) 针对电气传动控制系统中的直流电动机转速控制问题，应用“转速电流双闭环控制方案”解决了直流电动机的“稳速运行与最佳启动”的控制问题； 仿真实验结果证明了该方案的有效性。

(2) 将“双闭环控制策略”与“PID控制器的工程设计方法”应用于一阶直线倒立摆的控制系统设计中，仿真实验表明： 控制系统具有良好的动态性能与鲁棒性。

(3) 一阶并联旋转双倒立摆系统(PRDIP系统)是一个自不稳定的“非最小相位系统”，需要通过控制系统设计的方法来实现“闭环控制系统”的稳定； 当PRDIP系统的两根摆杆参数完全一致时候，系统是“不能控”的； 系统旋转部件的“摩擦系数”也会影响系统的可控性。

(4) 针对“摆长时变条件下的吊车重物防摆控制”问题，应用滑模变结构控制理论，在保证系统鲁棒性的前提下，有效地提高了重物定位的效率。

(5) MBD(基于模型的设计)方法是控制系统设计的基本方法，对于复杂的机电一体化控制系统设计，应用Matlab/Simscape“可视化模型”实现的系统仿真实验具有“仿真精度高(建模过程AI自动化)、动态响应直观(运动姿态VR可视化)”的特点。

(6) 自平衡式两轮电动车是一个典型的“机电一体化系统”，对其进行的系统建模与控制系统设计，有助于我们研究更有效的系统控制问题。

(7) 灵长类仿生机器人是一个典型的“欠驱动机械系统”，对其深入的探究有助于我们提高分析问题与解决问题的能力。

（8） 本章所涉及的电力电子、运动控制等内容，参考了相关领域的经典著作［3843］，建议读者尽可能阅读这些原著，以便于对内容的理解与深入。

需要说明的是，本章各节内容中的具体技术路线与结果，不一定是最佳(或者说是正确)的。但是，其工程背景与所归结的问题应该说都是值得我们深入探讨与研究的。因此，笔者希望通过本章的内容，能够起到一个抛砖引玉的作用，为读者提供一个畅想的空间，以达到“分析与解决复杂工程问题”能力培养的目的。


习题

31如图3147所示一带有库仑摩擦的二阶随动系统，试优化设计K1参数，并分析非线性环节对系统动态性能的影响。



图3147题31图


32试分析图3148所示系统中死区非线性对系统动态性能的影响。



图3148题32图


33如图3149所示计算机控制系统，试设计一最小拍控制器D(z)，并用仿真的方法分析最小拍控制器对系统输入信号和对象参数变化的适应性。



图3149题33图


34为使图3150所示系统不产生自激振荡，试分析a、b取值。



图3150题34图



35已知某一地区在有病菌传染下的描述三种类型人数变化的动态模型为


X·1=-αX1X2X1(0)=620
X·2=αX1X2-βX2X2(0)=10
X·3=βX2X3(0)=70


式中，X1表示可能传染的人数； X2表示已经得病的人数； X3表示已经治愈的人数； α=0.001;β=0.072。
试用仿真方法求未来20年内三种人人数的动态变化情况。

36如图3151所示旋转倒立摆系统，试讨论如下问题： 



图3151题46旋转倒立摆系统原理图


(1)　能否通过控制直流力矩电动机的输出转矩来实现“自由旋转摆”的垂直倒立(即θ1=θ2≈0)?

(2) 能否实现“自由旋转摆”在平面的适当移动(即θ1≈c(常数)，θ2≈0)?

(3) 试给出你的具体实现方案。

37在如图3152所示的某高精度齿轮测量机系统中，为保证主轴系统的回转精度，采用了“过盈量轴承”的装配技术，从而使主轴控制电动机在低速下的机械特性呈现图(b)所示的情况(带有随机因素影响的非线性库仑摩擦特性)。若电动机选用低速大惯量永磁力矩电动机(基本参数： 额定电压36V，额定转速30r/min、瞬时最大转矩50N·m)，回转角度检测选用分辨率为0.25″(角度)的高精度圆光栅，试分析回答如下问题： (1)若系统要求控制电动机达到(1r/10h~1r/min)的调速范围，试设计电控系统的主回路(即电动机的驱动器)、数字式给定电路以及数字控制器； (2)若测量工艺要求伺服系统定位精度达到±1个脉冲(0.25″)，定位步长为0.1°，试利用数字仿真技术设计该伺服系统，使得其动态性能具有“快速无超调”的特点。



图3152题37高精度齿轮测量机主轴系统结构图


38在大型立体化仓库中经常采用如图3153所示的搬运机器人，由于H值较高，水平方向位置检测码盘又装在顶部，所以在定位过程中往往会产生“位置抖动”(用θ来描述)现象，从而影响定位的精度，降低了定位过程的效率。试针对这一问题研究如何建立最佳控制规律u(t)(运行速度)，以使机器人从A点运行定位到B点所用的时间最短，同时又满足一定的定位精度要求。



问题探究
立体仓库
搬运机器
人控制问
题（机械手
末端防抖）



图3153题38立体仓库搬运机器人结构图


(提示： 本问题所讨论的问题可抽象为图3154所示的物理模型——具有弹性立杆的移动小车问题)



图3154具有弹性立杆的移动小车问题示意图



39如今，在没有深水港口的情况下，船载吊车广泛地用于大型集装箱运输船只的装卸以及军事舰艇的物料补给等方面(如图3155所示)。但是，由于吊运过程中随机海浪引起的负载摆动不但影响吊装工作的效率与准确性，而且还可能损坏货物，甚至引起人员伤亡。因此，负载防摆控制工作就显得尤为重要。




问题探究
船用吊车的
重物防摆控
制问题（船
舶的海浪摇
摆状况）




图3155题39船载吊车工作示意图



试根据这一问题，给出一种简单易行的负载消摆控制方案，使其可有效地抑制在下放货物时海浪对船体扰动所引起的货物摆动，并利用仿真技术加以验证。


310三相四开关PWM整流器作为三相六开关PWM整流器的容错拓扑，近年来成为学术研究热点问题，其拓扑结构如图3156所示。图中，ea,eb,ec为三相交流电压；  ia,ib,ic为交流侧电流；  Sa，Sb，及S′a,S′b为ab两相上下桥臂的功率开关管L；  L为网侧滤波电感；  R为交流侧等效电阻； C1,C2为直流侧两个滤波电容； RL为负载电阻； uo为直流侧电压。

试以系统高功率因数运行为控制目标，设计三相四开关PWM整流器的控制系统，并比较其与三相六开关PWM整流器控制系统设计的异同。



图3156三相四开关PWM整流器拓扑结构



参考文献

［1］阮毅，陈伯时.电力拖动自动控制系统——运动控制系统(第4版)［M］.北京: 机械工业出版社，2010.

［2］Felix Grasser.A Mobile，Inverted pendulum［J］.IEEE Transactions on industrial electronics,2002，49(1): 107114.

［3］屠运武，徐俊艳，张培仁，等．自平衡控制系统的建模与仿真［J］．系统仿真学报，2004，16(4)： 839841．

［4］阮晓钢，蔡建羡，李欣源，等．两轮自平衡两轮车的研究与设计［M］．北京： 科学出版社，2012： 110146．

［5］张晓华，张志军．自平衡式两轮电动车耦合控制研究［J］．控制工程，2013，20(1)： 2629．

［6］夏德黔，翁贻方．自动控制理论［M］．北京： 机械工业出版社，2004．

［7］M.J.Er，M.Zribi and K.L.Lee.Variable Structure Control of Overhead Crane［J］.Proceedings of the 1998 IEEE International Conference on Control  Applications Trieste，Ital 14 September,1998.

［8］向博,高丙团,张晓华.桥式吊车的滑模变结构控制［J］.控制工程，2006，6(23).

［9］高为炳，程勉.变结构控制的品质控制［J］.控制与决策，1989，4(4).

［10］向博,高丙团,张晓华,等.非连续系统的Simulink仿真方法研究［J］.系统仿真学报，2006，5(32).

［11］邱杰，原渭兰.数字计算机仿真中消除代数环问题的研究［J］.计算机仿真，2003，20(7).

［12］Kent H.Lundberg，James K.Roberge.Classical DualInvertedPendulum Control Proceedings of the IEEE CDC，2003：  43994404.

［13］Jianqiang Yia，Naoyoshi Yubazakib，Kaoru Hirotac.A new fuzzy controller for stabilization of paralleltype double inverted pendulum system［J］.Fuzzy Sets and Systems，2002，126： 105119.

［14］周国华，许建平.开关变换器调制与控制技术综述［J］.中国电机工程学报，2014，34(6): 815831.

［15］陆益民，张波，尹丽云.DC/DC变换器的切换仿射线性系统模型及控制［J］.中国电机工程学报，2008，28(15): 1622.

［16］卢伟国，赵乃宽，郎爽，等.Lyapunov模式控制Buck变换器［J］.电工技术学报，2014，29(10): 98105.

［17］Theunisse T，Chai J，Sanfelice R G，et al.Robust Global Stabilization of the DCDC Boost Converter via Hybrid Control［J］.IEEE Transactions on Circuits and Systems I: Regular Papers，2015，62(4): 10521061.

［18］Singh S，Fulwani D，Kumar V.Robust Slidingmode Control of DC/DC Boost Converter Feeding a Constant Power Load［J］.IET Power Electronics，2015，8(7): 12301237.

［19］Erickson R W，Maksimovic D.Fundamentals of Power Electronics［M］.Springer Science & Business Media，2001.

［20］徐德鸿.电力电子系统建模及控制［M］.北京： 机械工业出版社，2006.

［21］张卫平.开关变换器的建模与控制［M］.北京： 中国电力出版社，2005.

［22］Bing Z, Du X，Sun J.Control of Threephase PWM Rectifiers Using a Single DC Current Sensor［J］.IEEE Transactions on Power Electronics，2011，26(6): 18001808.

［23］Ge J, Zhao Z，Yuan L，et al.Direct Power Control Based on Natural Switching Surface for Threephase PWM Rectifiers［J］.IEEE Transactions on Power Electronics，2015，30(6): 29182922.

［24］Hemdani A，Dagbagi M，Naouar W M，et al.Indirect Sliding Mode Power Control for Three Phase Grid Connected Power Converter［J］.IET Power Electronics，2015，8(6): 977985.

［25］Aguilera R P，Quevedo D E.Predictive Control of Power Converters: Designs with Guaranteed Performance［J］.IEEE Transactions on Industrial Informatics，2015，11(1): 5363.

［26］Ketzer M B，Jacobina C B.Sensorless Control Technique for PWM Rectifiers with Voltage Disturbance Rejection and Adaptive Power Factor［J］.IEEE Transactions on Industrial Electronics，2015，62(2): 11401151.

［27］Lee Tzann Shin. Inputoutput Linearization and Zerodynamics Control of Threephase AC/DC Voltagesource Converters［J］.IEEE Transactions on Power Electronics，2003，18(1): 1122.

［28］Yin B，Oruganti R，Panda S K，et al.A Simple Singleinput—singleoutput (SISO) Model for a Threephase PWM Rectifier［J］.IEEE Transactions on Power Electronics，2009，24(3): 620631.

［29］张兴,张崇巍.PWM整流器及其控制(第2版)［M］.北京： 机械工业出版社，2012.

［30］杨德刚,赵良炳,刘润生.三相高功率因数整流器的建模及闭环控制［J］.电力电子技术,1999，33(5): 4951.

［31］朱永亮,马惠,张宗濂.三相高功率因数PWM整流器双闭环控制系统设计［J］.电力自动化设备,2006，26(11): 8791.

［32］Dannehl J, Wessels C，Fuchs F W.Limitations of Voltageoriented PI Current Control of Gridconnected PWM Rectifiers with Filters［J］.IEEE Transactions on Industrial Electronics，2009，56(2): 380388.

［33］H.Song, K.Nam.Dual Current Control Scheme for PWM Converter under Unbalanced Input Voltage Conditions［J］.IEEE Transactions on Industrial Electronics，1999，46(5): 953959.

［34］Zhang Y, Qu C.Model Predictive Direct Power Control of PWM Rectifier under Unbalanced Network Conditions［J］.IEEE Transactions on Industrial Electronics，2015，62(7): 40114022.

［35］Roiu D，Bojoi R I，Limongi L R，et al.New Stationary Frame Control Scheme for Threephase PWM Rectifiers under Unbalanced Voltage Dips Conditions［J］.IEEE Transactions on Industry Applications，2010，46(1): 268277.

［36］Valouch V, Bejvl M，Simek P，et al.Power Control of Grid Connected Converters under Unbalanced Voltage Conditions［J］.IEEE Transactions on Industrial Electronics，2015，62(7): 42414248.

［37］郭小强，李建，张学，等.电网电压畸变不平衡情况下三相PWM整流器无锁相环直流母线恒压控制策略［J］.中国电机工程学报，2015，35(8): 20022008.

［38］Bimal K Bose.Modern Power Electronics and AC Drives［M］.Prentice Hall，2001.

［39］王兆安，刘进军.电力电子技术(第5版)［M］.北京: 机械工业出版社，2009.

［40］(美)Richard C.Dorf，Robert H.Bishop.现代控制系统.第九版，北京： 科学出版社，2002.

［41］(法)Mohand Mokhtari，Michel Marie.MATLAB与Simulink工程应用.赵彦玲，等.北京： 电子工业出版社，2002.

［42］(日)细江繁幸.系统与控制.北京： 科学出版社，2001.

［43］(日)末松良一.机械控制入门.北京： 科学出版社，2000.

［44］Bingtuan G,Junyuan L,Jianguo Z.Modeling and stabilization control of an underactuated double rotating pendulums system： 2010 8th World Congress on Intelligent Control and Automation［C］.Jinan,2010： 159164.

［45］虞俊豪.一阶并联旋转双倒立摆系统的EFC/LQR双模态控制［D］.大连： 大连理工大学,2022.

［46］于振宝.一阶并联旋转双倒立摆系统全域运动控制［D］.大连： 大连理工大学，2023.

［47］哈尔滨工业大学理论力学教研室.理论力学（第9版）［M］.北京： 高等教育出版社,2023.

［48］王永刚.分析力学［M］.北京： 清华大学出版社,2019.