论文部分内容阅读
【摘要】本文结合“卓越工程师”的培养目标和本科生的学习特点，介绍“系统建模与仿真”课程的本科教学实践体会，借助于MATLAB软件用实例分析三种阶跃响应辨识方法，有助于学生深入理解建模原理的本质，培养学生的动手能力和创新能力，从而提高教学质量。
【关键词】本科；系统建模与仿真；教学实践；MATLAB
【Abstract】According to the requirements of outstanding engineers’ plan and learning characteristics of undergraduate students， the paper discusses some teaching methods of the course “System Modeling and Simulation”. As a teaching example， the paper introduces three kinds of step response identification methods and their realization by using MATLAB. The practice helps students understand in depth the essence of the modeling principles， develop practical and innovative ability， thus improves the teaching quality.
【Keywords】undergraduate；system modeling and simulation；teaching practice；MATLAB
1 引言
系统建模与仿真课程是许多高校自动化专业研究生的基础课程之一，也常被作为自动化高年级本科生的一门专业拓展课程[2]。该课程理论性强，内容丰富，同时与生产实践又紧密相连，应用范围很广。2010 年，教育部发布“卓越工程师培养计划”，旨在培养造就一大批实践创新能力强、适应社会发展需要的各类型高质量工程技术人才。在教学中如何针对本科生的特点，尤其是从培养“卓越工程师”的需求出发，把握好课程内容的重点和难点，培养学习的主观能动性、动手能力和创新能力，本文以教学中的一个具体实例对教学实践作了一些探讨。
2 课程特点及教学目的
系统建模与仿真是研究建立系统模型的一套理论和方法。内容涵盖建模、仿真和可信度评估三大方面。该课程的特点之一是理论性较强，涉及概率与统计，随机过程等多方面的数理知识，有些理论又比较抽象，是一门具有一定深度和难度的课程。该课程的第二个特点是内容丰富，涉及面广。国内外与该课程相关的课程主要有《系统辨识》、《系统建模与仿真》、《系统建模、辨识及仿真》、《系统建模与参数估计》等。目前，各学校在基于系统辨识的参数估计教学上，内容比较统一，以经典辨识方法、最小二乘法和极大似然估计法为主。但是在其他内容的教学上比较零散，每个学校侧重点均有不同[2]。有的以离散事件系统的建模与仿真为主，有的以生产物流系统作为主要对象，有的讲授基于Matlab/Simulink的控制系统的建模与仿真，有的则将系统辨识与自适应控制内容相结合。另外，该课程的第三个特点是具有很强的实践性，强调理论联系实际，必须通过实践教学环节才能使学生很好地理解和掌握系统建模与仿真的理论和方法。
东华大学针对信息学院的本科生在大四第一学期开设了该课程，课程的名称是《系统建模仿真与实践》。该课程首先要求学生对系统建模与仿真这一学科的体系架构、主要的研究内容和方向有个总体的了解和把握；其次，要求学生理解系统辨识的基本概念，掌握常用的系统辨识方法，淡化抽象的理论推导；最后，强调实践和应用，要求学生通过大量的教学案例和计算机仿真实验对基本原理和方法做到全面灵活的掌握。
3 教学实践举例
在教学过程中发现，教师在讲解基本原理和方法时，学生尚较容易理解。但是当让他们将原理变成算法在计算机上实现时，往往显得手足无措。一方面是因为对软件的熟悉程度不够，或者编程能力不强；另一方面则是因为对基本原理的理解还不够透彻。这就需要教师亲自示范，最好做到每讲解一个建模方法，就附加一个案例和相应的程序，借助仿真软件直观、便捷的特点帮助学生将基本原理和应用结合起来，提高学生的学习兴趣，并真正达到学以致用的目的。下面举例说明。
已知某系统阶跃响应的实验值如表1所示，试辨识得到系统的传递函数模型。
表1：某系统阶跃响应的实验值[3]
t y*（t）
6 0.86318
8 0.94516
10 0.97871
12 0.99189
14 0.99695
16 0.99886
18 0.99958
20 0.99984
解：首先根据实验数据画出系统的阶跃响应曲线，用MATLAB编写程序如下：
clear
clc
% 表1原始数据
t = [6 8 10 12 14 16 18 20]；y = [0.86318 0.94516 0.97871 0.99189 0.99695 0.99886 0.99958 0.99984]；
L=length（t）；
% 绘制实际系统阶跃响应曲线
figure（1）
plot（t，y）
title（'Step response of the real system'）
xlabel（'Time（sec）'） grid
图1：实际系统阶跃响应曲线
阶跃响应曲线如图1所示。根据该曲线的特点，选取二阶振荡系统作为系统的模型结构，如下所示：
G（s）=1T2s2+2ξTs+1 ，
（1）
其中，ξ>1，系统为过阻尼系统。由此，系统具有两个实特征根，写成零极点的形式为：
G（s）=ω1ω2（s+ω1）（s+ω2） ，
（2）
其中，
T=1ω1ω2
ξ=ω1ω2ω1ω2
对（2）式进行Laplace逆变换，得到系统的阶跃响应为：
y（t）=1-ω2ω1-ω2e-ω1t+ω1ω1-ω2e-ω2t
，（3）
令ω2=aω1，带入（3）式得
y（t）=1-aa-1e-ω2t+1a-1e-aω1t，
（4）
于是，
1-y（t）=aa-1e-ω1t（1-1ae-（a-1）ω1t），
（5）
方法一：线性拟合
将式（5）等式两边取自然对数，得到
ln（1-y（t））=lnaa-1-ω1t+ln（1-1ae-（a-1）ω1t），
（6）
当t→∞时，ln（1-1ae-（a-1）ω1t）→0，则（6）式简化为
ln（1-y（t））=lnaa-1-ω1t （t→∞），
（7）
（7）式意味着当t→∞时，ln（1-y（t））近似为一条直线，即
ln（1-y（t））=kt+b，
其中，
k=-ω1
b=lnaa-1，
（8）
所以，只要选取t较大时的实验数据ln（1-y（t））进行一次多项式拟合，得到渐近线的斜率k和截距b，就可求得ω1和a，进而求得系统的传递函数模型（1）。
将上述原理用MATLAB编写程序如下：
% 渐近线线性拟合
tm = t（L-3：L）；% t goes to infinity
ym = y（L-3：L）；
ym1 = log（1-ym）；
aa = polyfit（tm，ym1，1）；
k = aa（1）；
b = aa（2）；
% 拟合效果对比
figure（2）
plot（tm， k*tm+b，'r'，tm，ym1）
title（'Linear fit'）
grid
legend（'拟合直线'，'ln（1-y）'）
% 参数和系统模型求解？
w1 = -k；
a = exp（b）/（exp（b）-1）；
w2 = a*w1；
T = 1/sqrt（w1*w2）；
yita = （w1 + w2）/2 * T；
z = [ ]；
p =[-w1 ； -w2]；
k1 = w1 * w2；
sys = zpk（z，p，k1） % 系统辨识模型
% 绘制辨识模型的阶跃响应曲线
figure（3）
step（sys）
grid
运行上述程序，得到系统的传递函数模型为：
Zero/pole/gain：
0.36383
________________________________________
（s+0.4922）（s+0.7392）
绘制的曲线参见图2-3所示。
图2：线性拟合效果
图3：辨识模型的阶跃响应曲线
从仿真结果看，线性拟合非常吻合，辨识得到的传递函数模型接近真实系统。然而，在上述程序中，使用的核心函数是polyfit（），它是将t→∞的数据ln（1-y（t））在最小二乘原则下进行拟合，也就是说辨识得到的系统模型是以最接近ln（1-y（t））为前提的，而不是题目最初给出的实验数据y（t）。为此，可改用下面一种方法。
方法二：非线性拟合
观察（5）式发现，当t→∞时，1-1ae-（a-1）ω1t→1，则（5）式可化简为
1-y（t）=aa-1e-ω1t t→∞，
（9）
（9）式意味着当t→∞时，1-y（t）可以拟合成一条负指数的曲线，即
y（t）=1-ce-ω1t，
（10）
其中，
c=aa-1
将上述原理用MATLAB编写程序如下，代替上述程序中的"渐近线线性拟合"和"拟合效果对比"两部分，从而得到系统参数a和ω1：
% 非线性拟合
tm = t（L-3：L）；% t goes to infinity
ym = y（L-3：L）；
m0 = [1 1]；
m = lsqcurvefit（@chap3ex1fun，m0，t，y）；
a = m（1）/（m（1）-1）；
w1 = m（2）；
其中，子函数chap3ex1fun程序为：
function f=chap3ex1fun（m，t） f = 1-m（1）*exp（-m（2）*t）；
完整运行程序，得到系统的传递函数模型为：
Zero/pole/gain：
0.39271
________________________________________
（s+0.4624）（s+0.8493）
从仿真结果来看，方法二与方法一得到的系统传递函数模型略有不同，这是因为方法二采用的核心函数是lsqcurvefit（ ），它是在最小二乘原则下以最接近实验数据y（t）为前提的（参见eq. （10））。
方法一和二都是直接使用了系统辨识工具箱中的专有拟合函数，这些拟合函数都是基于最小二乘原则的。为了让学生更深入理解并能灵活使用最小二乘法，可以自行编写程序。具体方法如下。
方法三：最小二乘估计
沿用上述例子，在线性拟合中，每个观测数据点与拟合直线的偏差为：
e（ti）=yy（ti）-（kti+b） i=1，2，…，n
其中，yy（ti）=ln（1-y*（ti））可根据表1的实验数据得到。偏差的平方和为：
j（k，b）=∑ni=1e2（ti）
根据最小二乘原理，应选取合适的k和b，使得j（k，b）取极小值，即k和b应满足jk=0和jb=0。由此求得：
k=n∑n
i=1tiyy（ti）-∑n
i=1ti∑n
i=1yy（ti）n∑n
i=1t2
1-∑n
i=1ti∑n
i=1ti
b=∑n
i=1yy（ti）-k∑n
i=1tin
将上述原理用MATLAB编写程序如下：
% 最小二乘线性拟合
tm = t（L-3：L）； % t goes to infinity
ym = y（L-3：L）；
y1 = log（1-ym）；
n = 4； Lx = 0； Ly = 0； Lxx = 0； Lxy = 0； Lyy = 0；
for i = 1： 4
Lx = Lx + tm（i）；
Ly = Ly + y1（i）；
Lxx = Lxx + tm（i）*tm（i）；
Lxy = Lxy + tm（i）*y1（i）；
Lyy = Lyy + y1（i）*y1（i）；
end
k = （n*Lxy-Ly*Lx）/（n*Lxx-Lx*Lx）；
b = （Ly-k*Lx）/n；
完整运行程序，得到系统的传递函数模型与上述用polyfit函数拟合得到模型完全一致。
综上所述，在教学中采用了一个只有8个实验数据点的简单例子，分别运用三种方法辨识得到系统的传递函数模型，辨识原理清晰，编程实现直观易行。不仅可以使学生将抽象的辨识原理与应用结合起来，还能激发学生的学习兴趣和主观能动性，学会举一反三，从而增强分析问题和解决问题的能力。
4 结束语
系统建模与仿真是一门理论性强、内容丰富、应用面广的课程。要针对学生的特点，选取合适的教学内容，并在教学过程中注重理论联系实际，加强案例教学和实践环节。以上述教学实践为例，在突出教学重点的基础上，利用MATLAB仿真手段对几种建模方法进行综合比较，有助于提高学生的理论学习能力和工程实践能力，取得了不错的教学效果，符合“卓越工程师”的培养要求。
参考文献
[1] 刘娣， 朱松青， 林健， 应用型本科院校系统辨识课程教学改革探索[J]， 科技创新导报， No. 31， 2012： 174-175.
[2] 谢美华， 严奉霞， 系统科学专业《系统建模与参数估计》课程教学改革[J]， 高等教育研究学报， Vol. 33， No. 4， 2010： 107-109.
[3] 宋志安， 王文馨， 由阶跃响应曲线辨识传递函数的图解方法[J]， 山东科技大学学报（ 自然科学版） ， Vol. 22， No. 1， 2003： 61-63.
[作者简介]
陈亮（1976-），女（回），江苏常熟人，副教授，博士，主要研究方向为控制理论与应用；石红瑞（1968-），女，山西朔州人。副教授，博士，主要研究方向为控制理论与应用。
【关键词】本科；系统建模与仿真；教学实践；MATLAB
【Abstract】According to the requirements of outstanding engineers’ plan and learning characteristics of undergraduate students， the paper discusses some teaching methods of the course “System Modeling and Simulation”. As a teaching example， the paper introduces three kinds of step response identification methods and their realization by using MATLAB. The practice helps students understand in depth the essence of the modeling principles， develop practical and innovative ability， thus improves the teaching quality.
【Keywords】undergraduate；system modeling and simulation；teaching practice；MATLAB
1 引言
系统建模与仿真课程是许多高校自动化专业研究生的基础课程之一，也常被作为自动化高年级本科生的一门专业拓展课程[2]。该课程理论性强，内容丰富，同时与生产实践又紧密相连，应用范围很广。2010 年，教育部发布“卓越工程师培养计划”，旨在培养造就一大批实践创新能力强、适应社会发展需要的各类型高质量工程技术人才。在教学中如何针对本科生的特点，尤其是从培养“卓越工程师”的需求出发，把握好课程内容的重点和难点，培养学习的主观能动性、动手能力和创新能力，本文以教学中的一个具体实例对教学实践作了一些探讨。
2 课程特点及教学目的
系统建模与仿真是研究建立系统模型的一套理论和方法。内容涵盖建模、仿真和可信度评估三大方面。该课程的特点之一是理论性较强，涉及概率与统计，随机过程等多方面的数理知识，有些理论又比较抽象，是一门具有一定深度和难度的课程。该课程的第二个特点是内容丰富，涉及面广。国内外与该课程相关的课程主要有《系统辨识》、《系统建模与仿真》、《系统建模、辨识及仿真》、《系统建模与参数估计》等。目前，各学校在基于系统辨识的参数估计教学上，内容比较统一，以经典辨识方法、最小二乘法和极大似然估计法为主。但是在其他内容的教学上比较零散，每个学校侧重点均有不同[2]。有的以离散事件系统的建模与仿真为主，有的以生产物流系统作为主要对象，有的讲授基于Matlab/Simulink的控制系统的建模与仿真，有的则将系统辨识与自适应控制内容相结合。另外，该课程的第三个特点是具有很强的实践性，强调理论联系实际，必须通过实践教学环节才能使学生很好地理解和掌握系统建模与仿真的理论和方法。
东华大学针对信息学院的本科生在大四第一学期开设了该课程，课程的名称是《系统建模仿真与实践》。该课程首先要求学生对系统建模与仿真这一学科的体系架构、主要的研究内容和方向有个总体的了解和把握；其次，要求学生理解系统辨识的基本概念，掌握常用的系统辨识方法，淡化抽象的理论推导；最后，强调实践和应用，要求学生通过大量的教学案例和计算机仿真实验对基本原理和方法做到全面灵活的掌握。
3 教学实践举例
在教学过程中发现，教师在讲解基本原理和方法时，学生尚较容易理解。但是当让他们将原理变成算法在计算机上实现时，往往显得手足无措。一方面是因为对软件的熟悉程度不够，或者编程能力不强；另一方面则是因为对基本原理的理解还不够透彻。这就需要教师亲自示范，最好做到每讲解一个建模方法，就附加一个案例和相应的程序，借助仿真软件直观、便捷的特点帮助学生将基本原理和应用结合起来，提高学生的学习兴趣，并真正达到学以致用的目的。下面举例说明。
已知某系统阶跃响应的实验值如表1所示，试辨识得到系统的传递函数模型。
表1：某系统阶跃响应的实验值[3]
t y*（t）
6 0.86318
8 0.94516
10 0.97871
12 0.99189
14 0.99695
16 0.99886
18 0.99958
20 0.99984
解：首先根据实验数据画出系统的阶跃响应曲线，用MATLAB编写程序如下：
clear
clc
% 表1原始数据
t = [6 8 10 12 14 16 18 20]；y = [0.86318 0.94516 0.97871 0.99189 0.99695 0.99886 0.99958 0.99984]；
L=length（t）；
% 绘制实际系统阶跃响应曲线
figure（1）
plot（t，y）
title（'Step response of the real system'）
xlabel（'Time（sec）'） grid
图1：实际系统阶跃响应曲线
阶跃响应曲线如图1所示。根据该曲线的特点，选取二阶振荡系统作为系统的模型结构，如下所示：
G（s）=1T2s2+2ξTs+1 ，
（1）
其中，ξ>1，系统为过阻尼系统。由此，系统具有两个实特征根，写成零极点的形式为：
G（s）=ω1ω2（s+ω1）（s+ω2） ，
（2）
其中，
T=1ω1ω2
ξ=ω1ω2ω1ω2
对（2）式进行Laplace逆变换，得到系统的阶跃响应为：
y（t）=1-ω2ω1-ω2e-ω1t+ω1ω1-ω2e-ω2t
，（3）
令ω2=aω1，带入（3）式得
y（t）=1-aa-1e-ω2t+1a-1e-aω1t，
（4）
于是，
1-y（t）=aa-1e-ω1t（1-1ae-（a-1）ω1t），
（5）
方法一：线性拟合
将式（5）等式两边取自然对数，得到
ln（1-y（t））=lnaa-1-ω1t+ln（1-1ae-（a-1）ω1t），
（6）
当t→∞时，ln（1-1ae-（a-1）ω1t）→0，则（6）式简化为
ln（1-y（t））=lnaa-1-ω1t （t→∞），
（7）
（7）式意味着当t→∞时，ln（1-y（t））近似为一条直线，即
ln（1-y（t））=kt+b，
其中，
k=-ω1
b=lnaa-1，
（8）
所以，只要选取t较大时的实验数据ln（1-y（t））进行一次多项式拟合，得到渐近线的斜率k和截距b，就可求得ω1和a，进而求得系统的传递函数模型（1）。
将上述原理用MATLAB编写程序如下：
% 渐近线线性拟合
tm = t（L-3：L）；% t goes to infinity
ym = y（L-3：L）；
ym1 = log（1-ym）；
aa = polyfit（tm，ym1，1）；
k = aa（1）；
b = aa（2）；
% 拟合效果对比
figure（2）
plot（tm， k*tm+b，'r'，tm，ym1）
title（'Linear fit'）
grid
legend（'拟合直线'，'ln（1-y）'）
% 参数和系统模型求解？
w1 = -k；
a = exp（b）/（exp（b）-1）；
w2 = a*w1；
T = 1/sqrt（w1*w2）；
yita = （w1 + w2）/2 * T；
z = [ ]；
p =[-w1 ； -w2]；
k1 = w1 * w2；
sys = zpk（z，p，k1） % 系统辨识模型
% 绘制辨识模型的阶跃响应曲线
figure（3）
step（sys）
grid
运行上述程序，得到系统的传递函数模型为：
Zero/pole/gain：
0.36383
________________________________________
（s+0.4922）（s+0.7392）
绘制的曲线参见图2-3所示。
图2：线性拟合效果
图3：辨识模型的阶跃响应曲线
从仿真结果看，线性拟合非常吻合，辨识得到的传递函数模型接近真实系统。然而，在上述程序中，使用的核心函数是polyfit（），它是将t→∞的数据ln（1-y（t））在最小二乘原则下进行拟合，也就是说辨识得到的系统模型是以最接近ln（1-y（t））为前提的，而不是题目最初给出的实验数据y（t）。为此，可改用下面一种方法。
方法二：非线性拟合
观察（5）式发现，当t→∞时，1-1ae-（a-1）ω1t→1，则（5）式可化简为
1-y（t）=aa-1e-ω1t t→∞，
（9）
（9）式意味着当t→∞时，1-y（t）可以拟合成一条负指数的曲线，即
y（t）=1-ce-ω1t，
（10）
其中，
c=aa-1
将上述原理用MATLAB编写程序如下，代替上述程序中的"渐近线线性拟合"和"拟合效果对比"两部分，从而得到系统参数a和ω1：
% 非线性拟合
tm = t（L-3：L）；% t goes to infinity
ym = y（L-3：L）；
m0 = [1 1]；
m = lsqcurvefit（@chap3ex1fun，m0，t，y）；
a = m（1）/（m（1）-1）；
w1 = m（2）；
其中，子函数chap3ex1fun程序为：
function f=chap3ex1fun（m，t） f = 1-m（1）*exp（-m（2）*t）；
完整运行程序，得到系统的传递函数模型为：
Zero/pole/gain：
0.39271
________________________________________
（s+0.4624）（s+0.8493）
从仿真结果来看，方法二与方法一得到的系统传递函数模型略有不同，这是因为方法二采用的核心函数是lsqcurvefit（ ），它是在最小二乘原则下以最接近实验数据y（t）为前提的（参见eq. （10））。
方法一和二都是直接使用了系统辨识工具箱中的专有拟合函数，这些拟合函数都是基于最小二乘原则的。为了让学生更深入理解并能灵活使用最小二乘法，可以自行编写程序。具体方法如下。
方法三：最小二乘估计
沿用上述例子，在线性拟合中，每个观测数据点与拟合直线的偏差为：
e（ti）=yy（ti）-（kti+b） i=1，2，…，n
其中，yy（ti）=ln（1-y*（ti））可根据表1的实验数据得到。偏差的平方和为：
j（k，b）=∑ni=1e2（ti）
根据最小二乘原理，应选取合适的k和b，使得j（k，b）取极小值，即k和b应满足jk=0和jb=0。由此求得：
k=n∑n
i=1tiyy（ti）-∑n
i=1ti∑n
i=1yy（ti）n∑n
i=1t2
1-∑n
i=1ti∑n
i=1ti
b=∑n
i=1yy（ti）-k∑n
i=1tin
将上述原理用MATLAB编写程序如下：
% 最小二乘线性拟合
tm = t（L-3：L）； % t goes to infinity
ym = y（L-3：L）；
y1 = log（1-ym）；
n = 4； Lx = 0； Ly = 0； Lxx = 0； Lxy = 0； Lyy = 0；
for i = 1： 4
Lx = Lx + tm（i）；
Ly = Ly + y1（i）；
Lxx = Lxx + tm（i）*tm（i）；
Lxy = Lxy + tm（i）*y1（i）；
Lyy = Lyy + y1（i）*y1（i）；
end
k = （n*Lxy-Ly*Lx）/（n*Lxx-Lx*Lx）；
b = （Ly-k*Lx）/n；
完整运行程序，得到系统的传递函数模型与上述用polyfit函数拟合得到模型完全一致。
综上所述，在教学中采用了一个只有8个实验数据点的简单例子，分别运用三种方法辨识得到系统的传递函数模型，辨识原理清晰，编程实现直观易行。不仅可以使学生将抽象的辨识原理与应用结合起来，还能激发学生的学习兴趣和主观能动性，学会举一反三，从而增强分析问题和解决问题的能力。
4 结束语
系统建模与仿真是一门理论性强、内容丰富、应用面广的课程。要针对学生的特点，选取合适的教学内容，并在教学过程中注重理论联系实际，加强案例教学和实践环节。以上述教学实践为例，在突出教学重点的基础上，利用MATLAB仿真手段对几种建模方法进行综合比较，有助于提高学生的理论学习能力和工程实践能力，取得了不错的教学效果，符合“卓越工程师”的培养要求。
参考文献
[1] 刘娣， 朱松青， 林健， 应用型本科院校系统辨识课程教学改革探索[J]， 科技创新导报， No. 31， 2012： 174-175.
[2] 谢美华， 严奉霞， 系统科学专业《系统建模与参数估计》课程教学改革[J]， 高等教育研究学报， Vol. 33， No. 4， 2010： 107-109.
[3] 宋志安， 王文馨， 由阶跃响应曲线辨识传递函数的图解方法[J]， 山东科技大学学报（ 自然科学版） ， Vol. 22， No. 1， 2003： 61-63.
[作者简介]
陈亮（1976-），女（回），江苏常熟人，副教授，博士，主要研究方向为控制理论与应用；石红瑞（1968-），女，山西朔州人。副教授，博士，主要研究方向为控制理论与应用。