第5章实例1： 打印机驱动控制系统分析与设计


打印机（Printer）是计算机的重要输出设备之一，是由约翰·沃特和戴夫·唐纳德合作发明的，用于将计算机的运算结果或中间结果以人所能识别的数字、字母、符号和图像等，依照规定的格式打印在相关介质上的设备。打印分辨率、打印速度和噪声是衡量打印机质量的重要指标。打印机的分类很多，如按打印元件对介质是否有击打动作，分为击打式打印机与非击打式打印机； 按打印字符结构，分为全形字打印机和点阵字符打印机； 按一行字在纸上形成的方式，分为串式打印机与行式打印机； 按采用的实现技术，分为喷墨式、激光式、热敏式、静电式、磁式、发光二极管式等打印机类型。目前，打印机正向轻、薄、短、小、低功耗、高速度和智能化方向发展。
打印机驱动系统是负责打印速度的重要单元，本章以打印机皮带驱动系统为例，采用MATLAB软件实现打印机驱动系统现代控制方法设计。首先，介绍打印机皮带驱动系统的状态空间模型，再对建立的状态空间模型分析打印机皮带驱动系统的稳定性、能控性、能观性等性能。进一步，运用极点配置状态反馈控制方法，设计打印机皮带驱动系统的状态反馈控制器和跟踪控制器，提高打印机皮带驱动系统的动态性能和鲁棒性能。最后，考虑打印机皮带驱动系统难以直接测量状态变量，设计系统的状态观测器，估计该系统中难以测量的状态变量，减少实际应用中传感器的使用数量，节约打印机生产成本和维护成本。

5.1打印机皮带驱动系统模型
5.1.1驱动系统微分方程


计算机常用的低价位打印机都有皮带驱动器，用来驱动打印设备沿打印页面横向移动。图5.1.1给出了一个装有直流电机的打印机皮带驱动系统示意图，其中，光传感器用来测量打印设备的位置，皮带张力变化用来调整滑轮的转动。


图5.1.1打印机皮带驱动系统示意图



为建立如图5.1.1所示的打印机皮带驱动系统的微分方程模型，设定打印机皮带驱动系统的皮带弹性系数为k，滑轮半径为r，电机轴转角为θ，右滑轮转角为θp，打印设备质量为m。打印设备在时间t的位置值为y(t)，可以通过光传感器测量位置变量y。光传感器的输出电压v1=k1y，控制器的输出电压为v2，v2与电机激励磁场相关。此外，控制器输出电压v2是传感器输出电压v1的函数，通常满足如下线性关系： 


v2=-k2dv1dt+k3v1(5.1.1)


其中，k2为微分系数，k3为比例系数。具体参数和受力分析如图5.1.2所示，其中，T1和T2分别为皮带上下侧的张力，θ 和θp分别为电机和滑轮的转角。令电机和滑轮的总转动惯量为J=J电机+J滑轮。通常采用中等功率的直流电机，如选定电机功率为典型的18马力（1马力=735.499W）情况下，有J=0.01kg·m2。进一步，假设磁场电感忽略不计，磁场电阻R=2Ω，电机常数Km=2N·m/A，而电机和滑轮的摩擦系数b=0.25N·ms/rad，滑轮半径r=0.15m。


图5.1.2打印机皮带驱动系统模型简图



现在列写打印机皮带驱动系统的运动方程。注意到y=rθp，于是皮带张力T1可表示为


T1=k(rθ-rθp)=k(rθ-y)(5.1.2)


张力T2为T2=k(y-rθ)，作用在质量m上的净张力为


T1-T2=md2ydt2(5.1.3)


且T1-T2=k(rθ-y)-k(y-rθ)=2k(rθ-y)=2kx1，其中第一个状态变量为x1=(rθ-y)。令第二个状态变量为x2=dy/dt，并利用式（5.1.2）和式（5.1.3）可得


dx2dt=2kmx1(5.1.4)


当选择第三个状态变量为x3=dθ/dt时，x1的一阶导数为


dx1dt=rdθdt-dydt=rx3-x2(5.1.5)


接下来，寻求描述电机旋转运动的微分方程。当L=0时，有电机磁场电流i=v2/R和电机转矩Tm=Kmi。因此


Tm=KmRv2(5.1.6)


而电机输出转矩提供了驱动皮带的转矩与扰动或者不希望的负载转矩之和，因此Tm=T+Td。转矩T驱动电机轴带动滑轮运动，满足


T=Jd2θdt2+bdθdt+r(T1-T2)(5.1.7)


因此


dx3dt=d2θdt2(5.1.8)


所以有


dx3dt=-Kmk1k2JRx2-bJx3-2krJx1-TdJ(5.1.9)


5.1.2驱动系统状态空间模型
式（5.1.4）、式（5.1.5）、式（5.1.9）是描述系统的3个一阶微分方程，其矩阵形式为


x·= 0-1 r

2km0 0

-2krJ-Kmk1k2JR-bJx+0
0
-1JTd(5.1.10)


假设某个打印机皮带驱动系统的结构参数取值如下： 质量m=0.2kg，光传感器k1=1V/m，半径r=0.15m，电感L=0，摩擦系数b=0.25N·ms/rad，电阻R=2Ω，常数Km=2N·m/A，惯量J=0.01kg·m2，弹性系数k=20和k2=0.1，且选择使用k2=0.1和k3=0（速度反馈）。把以上参数代入打印机皮带驱动系统的状态方程，可得


x·=0-10.15

20000

-600-10-25x+0

0

-100u(5.1.11)


定义系统输出变量为


y=-x3


则连续时间驱动系统的状态空间模型为


x·=0-10.15

20000

-600-10-25x+0

0

-100u

y=00-1x(5.1.12)



5.1.3驱动系统状态空间模型的离散化
以T=0.05s为采样周期，采用MATLAB软件中的c2d函数，将打印机皮带驱动系统状态空间模型离散化。在MATLAB命令窗口中输入如下指令： 

>> A=［0 -1 0.15;200 0 0;-600 -10 -25］; % 系统的状态矩阵

>> B=［0;0;-100］;% 系统的输入矩阵

>> ［G,H］=c2d(A,B,0.05) % 调用c2d函数


运行结果如下： 

G =

 0.6871 -0.04580.0037

 8.9114 0.75940.0241

-16.2917 0.23720.2359

H =

-0.0120

-0.0452

-2.7492


因此，所求的离散化状态空间模型为

x(k+1)=0.6871-0.04580.0037
8.91140.75940.0241
-16.29170.23720.2359x(k)+-0.0120
-0.0452
-2.7492u(k)
y(k）=00-1x(k)(5.1.13)



若取采样周期为T=0.01s，则类似可得相应的离散状态空间模型为


x(k+1)=0.9858-0.010.0013
1.99050.990.0014
-5.3739-0.06050.775x(k)+-0.0007
-0.0005
-0.8835u(k)
y(k)=00-1x(k)(5.1.14)


从以上两个离散状态空间模型可以看出，不同的采样周期所得到的离散模型是不同的。这也验证了离散状态空间模型依赖于所选取的采样周期的结论。

5.1.4驱动系统状态空间模型转换为传递函数模型
采用MATLAB软件中的ss2tf函数，将打印机皮带驱动系统状态空间模型转换为传递函数模型。在MATLAB命令窗口中输入如下指令： 

>> A=［0 -1 0.15;200 0 0;-600 -10 -25］;% 系统的状态矩阵

>> B=［0;0;-100］;% 系统的输入矩阵

>> C=［0 0 -1］;% 系统的输出矩阵

>> D=0; % 系统的直接转移矩阵

>> ［num,den］=ss2tf(A,B,C,D) % 调用ss2tf函数


运行结果如下： 

num =

1.0e+04 *

00.0100 02.0000

den =

1.0e+03 *

 0.00100.02500.29005.3000 


因此，所求系统的传递函数为


G(s)=100s2+20000s3+25s2+290s+5300(5.1.15)



5.2打印机皮带驱动系统性能分析
5.2.1驱动系统运动响应分析
1.  驱动系统单位脉冲响应分析



为了获得打印机皮带驱动系统的单位脉冲响应，在MATLAB软件中调用impulse函数，编写如下M文件（Example521.m）： 

% 系统的状态矩阵

A=［0 -1 0.15;200 0 0;-600 -10 -25］;

% 系统的输入矩阵

B=［0;0;-100］;

% 系统的输出矩阵

C=［0 0 -1］; 

% 系统的直接转移矩阵 

D=［0］; 

% 调用impulse函数

impulse(A,B,C,D) 

ylabel('y')

xlabel('时间/s')

grid on 


运行Example521.m文件，结果如图5.2.1所示。


图5.2.1系统的单位脉冲响应曲线




2.  驱动系统初始状态响应分析
为了获得打印机皮带驱动系统的初始状态响应，在MATLAB软件中调用initial函数，编写如下M文件（Example522.m）： 

% 系统的状态矩阵

A=［0 -1 0.15;200 0 0;-600 -10 -25］;

% 系统的输入矩阵

B=［0;0;-100］;

% 系统的输出矩阵

C=［0 0 -1］;

% 系统的直接转移矩阵 

D=［0］;

% 系统的初始状态

x0=［2,1,1］; 

% 调用initial函数

［y,x］=initial(A,B,C,D,x0); 

t=0:0.05:55.20;

%创建子图

subplot(3,1,1)

% 调用plot画图函数 

plot(t,x(:,1)) 

% 给x轴加标注 

xlabel('时间/s') 

% 给y轴加标注 

ylabel('x_1')

subplot(3,1,2)

plot(t,x(:,2))

xlabel('时间/s')

ylabel('x_2')

subplot(3,1,3)

plot(t,x(:,3))

xlabel('time/s')

ylabel('x_3')


运行Example522.m文件，结果如图5.2.2所示。


图5.2.2系统的初始状态响应曲线



5.2.2驱动系统能控性和能观性分析
（1） 系统状态能控的充分必要条件为


rankΓcA,B=rankBAB…An-1B=n(5.2.1)


即状态能控性判别矩阵Γc［A,B］满秩意味着系统是能控的，系统能控意味着系统出现偏差时，总可以通过控制使得偏差为零。
在MATLAB软件命令行窗口输入如下指令： 

% 系统的状态矩阵

>> A = ［0 -1 0.15; 200 0 0;-600 -10 -25］;

% 系统的输入矩阵

>> B=［0;0;-100］;

% 系统的输出矩阵

>> C=［0 0 -1］;

% 系统的直接转移矩阵

>> D=0;

% 获取能控性判别矩阵的秩

>> n=rank(ctrb(A,B)) 


运行结果如下： 

n =

 3


因此可得系统的状态能控性判别矩阵为


rank(ΓcA,B)=rank0-15375

00-3000

-1002500-53500=3(5.2.2)


可以看出，该状态能控性判别矩阵是满秩的，所以系统是能控的。
（2）  系统状态能观的充分必要条件为


rankΓoA,C=rankC

CA

︙

CAn-1=n(5.2.3)


在MATLAB软件命令行窗口输入如下指令： 

% 系统的状态矩阵

>> A = ［0 -1 0.15; 200 0 0;-600 -10 -25］;

% 系统的输入矩阵

>> B=［0;0;-100］;

% 系统的输出矩阵

>> C=［0 0 -1］;

% 系统的直接转移矩阵

>> D=0;

% 获取能控性判别矩阵的秩

>> n=rank(obsv(A,C)) 


运行结果如下： 

n =

 3


因此可得系统的状态能观性判别矩阵为


rank(ΓoA,C)=rank1.000000
0-1.00000.1500
-1002500-53500=3(5.2.4)


可以看出，该状态能观性判别矩阵是列满秩的（即矩阵的秩等于列数），所以系统是状态能观的，即可以通过系统的输出来观测系统中未知的状态，为观测器的设计提供了理论基础。

5.2.3驱动系统稳定性分析
对于该系统，采用Lyapunov方程处理方法，即线性时不变系统在平衡点处渐近稳定的充分必要条件是对任意给定的对称正定矩阵Q，存在一个对称正定矩阵P，使得方程


ATP+PA=-Q(5.2.5)


成立，那么该系统是渐近稳定的。
取Q为三阶单位矩阵，在MATLAB命令窗口中输入如下指令： 

% 系统的状态矩阵

>> A=［0 -1 0.15;200 0 0;-600 -10 -25］; 

% 对称正定的矩阵，这里取为三阶单位矩阵

>> Q=eye(3); 

% 调用lyap函数 

>>P=lyap(A',Q)

% 调用eig函数求特征值 

>> val=eig(P)


运行结果如下： 

P =

137.87370.63170.2114

0.63170.6604 -0.0132

0.2114 -0.01320.0213

val =

0.0206

0.6578

137.8769


由运算结果可得，矩阵P的特征值都是正的，所以矩阵P正定，所以系统在原点处的平衡状态是渐近稳定的。

5.3打印机皮带驱动系统控制器设计
5.3.1驱动系统极点配置状态反馈控制器设计


运用现代控制理论中的状态反馈和极点配置方法，对打印机皮带驱动系统设计状态反馈控制器。因为驱动系统是状态能控的，因此该系统可以通过状态反馈任意配置闭环极点。所设定的闭环极点可根据不同的性能指标进行修改。设计由式（5.1.10）描述的驱动系统的状态反馈控制器为u=-Kx，其中K的值由所配置的极点决定。
设计者希望通过引入状态反馈控制器使输出超调量σ≤5%，峰值时间tp≤0.5s。对此，设系统有一对主导极点，并且有


λ1,2=-ξωn±jωn1-ξ2(5.3.1)


其中，ξ和ωn是二阶系统的阻尼比和无阻尼自振频率。利用二阶系统的阻尼比ξ、无阻尼自振频率ωn等参数和超调量σ、峰值时间tp的关系，由


σ=exp-ξπ1-ξ2≤5%,tp=πωn1-ξ2≤0.5(5.3.2)


可得


ξ≥0.707,ωn≥9 


通过式（5.3.2）求得闭环系统的主导极点，而第3个极点选取在远离主导极点的左半开复平面中。因此，选取第3个极点为


pole=-7.07+7.07j-7.07-7.07j-100(5.3.3)


采用爱克曼公式法可以得到要设计的状态反馈控制器为


u=-81.59804.3770-0.8914x(5.3.4)


从而得到校正后的闭环系统为


x·=0-10.15
20000
-8759.8427.7-114.1x+0
0
-100u
y=00-1x(5.3.5)


在MATLAB软件中，编写如下M文件（Example531.m）： 

% 系统状态矩阵

A=［0 -1 0.15;200 0 0;-600 -10 -25］; 

% 系统输入矩阵

B=［0;0;-100］; 

% 系统输出矩阵 

C=［0 0 -1］;

% 系统直接转移矩阵 

D=［0］;

%系统的极点配置

J=［-100 -7.07+7.07*j -7.07-7.07*j］;

%调用acker函数求取状态反馈增益矩阵

K=acker(A,B,J); 

t=0:0.01:4;

% 调用step函数，求配置极点后的系统阶跃响应

［y,x1,t］=step(A-B*K,B,C,D,1,t); 

% 调用step函数，求原系统阶跃响应

［yy,x2,t］=step(A,B,C,D,1,t);

subplot(2,1,1)

plot(t,yy)

xlabel('时间/s')

ylabel('yy')

grid on

subplot(2,1,2);

plot(t,y)

grid on

xlabel('时间/s')

ylabel('y')


运行Example531.m文件，结果如图5.3.1所示。


图5.3.1开环系统与闭环系统的阶跃响应曲线




从图5.3.1可以看出，原系统阶跃响应有较为明显的振荡情况，从图5.3.1（a）可以大致判断其调节时间约为2.3s(Δ=5)，超调量为15%。从控制对象来考虑，打印机打印头的控制容错率非常低，因此无论从调节时间还是超调量来看，原系统动态性能都非常差。因此需要配置极点使闭环系统具有更好的过渡过程性能，在根据需求配置极点后，得到图5.3.1（b）的响应曲线。可以明显看到，配置极点后系统超调量非常小，调节时间不到0.5s，系统动态性能显著提升，但此时的稳态值不再是原系统阶跃响应稳态值3.774，出现了稳态误差。
5.3.2驱动系统跟踪控制器设计
由原打印机驱动系统的阶跃响应和极点配置后闭环系统的阶跃响应曲线图5.3.1可知，校正后的系统存在稳态误差，即极点配置方法可能会使一个原来没有稳态误差的系统产生稳态误差。另一方面，实际驱动系统不可避免地存在外部扰动，比如外力碰撞打印机等。对于其中的确定性扰动，扰动的存在使得打印机驱动系统在稳态时不能很好地跟踪参考输入，从而产生输出稳态误差。因此，需要对打印机驱动系统设计能够实现无静差跟踪阶跃参考输入信号的渐近跟踪调节器，即跟踪控制器。
首先，建立原打印机驱动系统的增广系统


x·
q·=A0
C0x
q+B
0u

y=C0x
q(5.3.6)


为了减小增广系统所增加的动态环节对原打印机驱动闭环系统性能的影响，要将增加的期望极点配置在打印机驱动闭环极点的左边，则再选择一个极点为-9。
在MATLAB软件中，编写如下M文件（Example532.m）： 

% 增广系统极点配置

A=［0 -1 0.15;200 0 0;-600 -10 -25］;B=［0;0;-100］;C=［0 0 -1］;D=［0］;

% 增广系统的状态矩阵

AA=［A zeros(3,1); C 0］;

% 增广系统的输入矩阵

BB=［B;0］;

% 增广系统极点配置

JJ=［-100 -7.07+7.07*1i -7.07-7.07*1i -9］;

% 调用acker函数求取状态反馈增益矩阵

KK=acker(AA,BB,JJ);

% 闭环系统参数矩阵

K1=［KK(1) KK(2) KK(3)］;K2=KK(4);

% 增广闭环系统状态矩阵

Ac=［A-B*K1 -B*K2;C 0］;

% 增广闭环系统输入矩阵

Bc=［0;0;0;-3.774］; 

% 增广闭环系统输出矩阵

Cc=［C 0］; 

% 增广闭环系统直接转移矩阵 

Dc=0;

% 绘制闭环输出响应曲线

t=0:0.01:4;

［yyy,x3,t］=step(Ac,Bc,Cc,Dc,1,t);

plot(t,yyy)

grid

xlabel('时间/s')

ylabel('y')


运行Example532.m文件，结果如图5.3.2所示。


图5.3.2加入跟踪控制器后闭环系统的阶跃响应曲线




在程序代码中，将跟踪的参考输出量设定为原系统阶跃响应稳态值3.774，这样可以方便与原系统阶跃响应曲线进行比较。事实上，在实际系统中可以任意调整参考输出量来控制最后输出稳态值。从图5.3.2所示曲线可以看出，改进后的系统完美继承了原来由配置极点方法所设计闭环系统的动态性能，而且消除了稳态误差，从而进一步改善了系统性能。值得一提的是，通过设计跟踪控制器，原先由极点配置控制器得到的阶跃响应曲线中无法消除的小尖峰也削除了，由此可得跟踪控制器增加的新极点对系统稳定性也有一定帮助。

5.3.3驱动系统状态观测器设计
对打印机皮带驱动系统模型进行分析，可以发现系统内部的状态信号不完全是可以从外部直接测量得到的，或者某些信号非常不容易测量，如角加速度难以直接通过传感器测量，而且引入过多的传感器也会在一定程度上降低控制系统的可靠性。所以直接进行状态反馈的控制是不合适的。
从系统的性能分析中可知该系统是状态能观的，那么它所有的状态信息都必定在输出中得到反映，因此可以设计一个状态观测器，用系统的外部输入输出信息来确定其状态量，避免直接测量打印机皮带驱动系统的状态变量。
考虑由式（5.1.10）描述的打印机皮带驱动系统，设计龙伯格观测器为


x=Ax～+Bu+L(y-Cx～)=(A-LC)x～+Bu+Ly(5.3.7) 


其中，误差e=x-x～的动态变化为


e·=(A-LC)e(5.3.8) 


在建立观测器后，可以进一步通过仿真得出误差曲线来检验观测器的效果。
在MATLAB命令窗口中输入如下指令： 

% 系统的状态矩阵

>> A=［0 -1 0.15;200 0 0;-600 -10 -25］;

% 系统的输出矩阵 

>> C=［0 0 -1］;

% V为观测器极点

>> V=［-200 -14.14+j*14.14 -14.14-j*14.14］; 

% L为观测器增益矩阵

>> L=(acker(A',C',V))';


得到观测器的增益矩阵为

L =

9.9993

-23.3688

-203.2800



且A-LC得

AA =

0-1.000010.1493

200.00000-23.3688

-600.0000-10.0000-228.2800


相应的观测器为


x=0-110.1493
2000-23.3688
-600-10-228.2800x～+0
0
-100u+9.9993
-23.3688
-203.2800y(5.3.9) 


状态估计的误差动态方程为


e·=(A-LC)e=0-110.1493
2000-23.3688
-600-10-228.2800e(5.3.10)


下面进一步通过仿真来检验观测器的效果。取初始误差为


e(0)=121.5T



在MATLAB软件中，编写如下M文件（Example534.m）： 

% AA为误差系统状态矩阵

AA=［0 -1 10.1493;200 0 -23.3688;-600 -10 -228.28］;

% BB为误差系统输入矩阵

BB=［0;0;0］;

C=［0 0 -1］;

D=0;

e0=［1;2;1.5］; % 初始误差向量

t=0:0.01:4;

sys=ss(AA,BB,C,D); % 误差的状态空间模型

% 获得误差系统的初始状态响应

［y,t,e］=initial(sys,e0,t);

% 画出状态变量1的误差曲线

subplot(3,1,1),plot(t,e(:,1)) 

grid

ylabel('e1')

xlabel('时间/s')

subplot(3,1,2),plot(t,e(:,2)) % 画出状态变量2的误差曲线

grid

ylabel('e2')

xlabel('时间/s')

subplot(3,1,3)

plot(t,e(:,3)) % 画出状态变量3的误差曲线

grid

ylabel('e3')

xlabel('时间/s') 


运行Example534.m文件，可得到状态估计的误差曲线如图5.3.3所示。


图5.3.3状态估计的误差曲线




从图5.3.3中可以看出，尽管系统的真实状态和观测器状态的初值有误差，但是随着时间的推移，它们之间的误差将衰减到零。由此可以得出，该观测器符合期望要求。采用这个观测器时，只需要测量输出x3即角速度就可以估计系统的所有状态信息。