第一部分控制系统仿真基础知识及 Matlab 在控制系统仿真中...

目录第一部分控制系统仿真基础知识及 Matlab 在控制系统仿真中的应用.1 实验一数学模型的转换.1 1.1. Matlab 符号数学介绍 1 1.2. 数学模型之间的转换.8 1.2.1 传递函数转换为可控标准型.8 1.2.2 传递函数转换为可观标准型.12 1.2.3 输入输出同阶次的处理方法.13 1.2.4 状态初值的通用求解方法.14 实验内容:14 实验二连续系统离散化.16 2.1 Matlab 基本函数 16 2.2 欧拉法.17 2.3 龙格-库塔法.20 2.4 亚当斯-巴什福思(Adams-Bashforth)22 2.5 预估-校正法.22 2.6 状态转移法离散化.23 2.7 连续系统传递函数到离散系统模型的转化.23 实验内容:26 实验三连续控制系统仿真.27 3.1 面向闭环系统传递函数的数字仿真.27 3.2 面向系统结构的数字仿真.27 3.2.1 与信息流无关的仿真方法.27 3.2.2 与信息流相关的仿真方法.29 3.3 非线性环节仿真.30 3.4 数字控制器仿真.33 实验内容:33 第二部分控制系统计算机辅助设计.35 实验四计算机辅助建立系统动态模型.35 4.1 经典时域辨识法.35 4.1.1 方波响应与阶跃响应之间的相互转换.35 4.1.2 利用阶跃响应数据获得对象的低阶传递函数模型.37 4.1.3 从响应曲线求对象的微分方程模型.43 4.2 频域辨识方法.45 4.2.1 由系统脉冲过渡函数计算频率特性.45 1 4.2.2 Levy 法拟合系统的传递函数.46 4.3 相关分析辨识系统脉冲响应函数.52 4.3.1 相关函数的求取.52 4.3.2 用伪随机信号辨识脉冲响应函数.53 实验内容:55 实验五控制系统计算机辅助分析与设计的频域方法.56 5.1 传递函数的频率特性计算.56 5.2 串联校正装置传递函数的计算.58 实验内容:59 实验六 PID 控制器的计算机仿真与辅助设计 60 6.1 连续系统 PID 控制器仿真及设计 60 6.2 数字 PID 控制器仿真及设计 64 实验内容:65 实验七基于时域的线性控制系统计算机辅助设计.66 6.1 极点配置.66 6.2 具有线性二次型性能指标的最优控制器设计.70 实验内容:71 综合实验.73 综合实验示例.73 实验内容:82 1 第一部分控制系统仿真基础知识及 Matlab 在控制系统仿真中的应用实验一数学模型的转换实验目的: 熟悉 Matlab 中的符号运算,掌握传递函数模型与状态空间模型之间的转换方法, 以及状态初值的求解方法. 实验方法: 1.1. Matlab 符号数学介绍 1) 创建符号变量和表达式 Matlab 符号数学工具箱提供了两个基本函数, 用来创建符号变量和表达式, 分别是 sym 和syms x=sym('x') 创建一个符号变量 x syms a b c 一次创建多个符号变量例: 2) 符号表达式的简化 ? 因式分解符号表达式的因式分解函数为 factor,调用格式: factor(S) 此函数因式分解符号表达式 S 的各个元素. 例: ? 符号表达式的展开符号表达式的展开函数为 expand,调用格式: expand(S) 此函数因式展开符号表达式 S >> syms x >> f=factor(x^9-1) f = (x-1)*(x^2+x+1)*(x^6+x^3+1) >> a=sym('a') a = a >> b=sym('hello') b = hello 2 例: ? 符号表达式的同类项合并 Matlab 的符号数学工具箱提供了两个化简函数,分别为 simple 和simplify,它们的调用格式: simplify(S) 化简函数,可用于化简各种表达式 [r,how]=simple(S) 寻找符号 S 的最简型, r 为返回的最简化形式, how 为化简过程所使用的主要方法. 例: ? 符号表达式的分式通分符号表达式的分式通分函数为 numden,其调用函数: [n,d]=numden(S) 将符号表达式 S 转换为分子和分母都是整系数的最佳多项式例:对表达式 x y y x f + = >> syms x y >> f=(x+1)^5; >> expand(f) ans = x^5+5*x^4+10*x^3+10*x^2+5*x+1 >> f=sin(x+y); >> expand(f) ans = sin(x)*cos(y)+cos(x)*sin(y) >> syms x >> f=sin(x)^2+cos(x)^2; >> simplify(f) ans = 1 >> f=exp(x*i)+exp(-x*i) f = exp(i*x)+exp(-i*x) >> simplify(f) ans = 2*cos(x) 3 ? 符号表达式的替换函数 subs 是用指定符号替换符号表达式中的某一特定符号,调用格式为: R=subs(S,old,new)用新的符号变量 new 替换原来符号表达式 S 中的变量 old.当变量 new 是数值形式时,显示的结果虽然是数值,但事实上是符号变量. 例: ? 符号表达式到多项式系数 SYM2POLY 符号表达式到多项式系数向量 SYM2POLY(P) 返回包含符号多项式 P 的系数的行向量 POLY2SYM 多项式系数向量到符号表达式 POLY2SYM(C) 返回系数为 C 的符号多项式(缺省符号为 x). Example: sym2poly(x^3 - 2*x - 5) returns [1 0 -2 -5]. COEFFS 返回多元(多变量)多项式系数. C = COEFFS(P) 返回多元(多变量)多项式 P 的系数 3) 符号微积分符号极限 limit(F,x,a) 计算符号表达式 F 在x→a 条件下的极限 limit(F,a) 计算符号表达式 F 中由默认自变量趋向于 a 条件下的极限 limit(F) 计算符号表达式 F 中由默认自变量趋向于 0 的极限 limit(F,x,a,'right'),limit(F,x,a,'left') 计算符号表达式 F 在x→a 条件下的右/左极限例:分别计算表达式 x x x ) sin( lim 0 → , x x 1 lim 0+ → , x x x ) 1 1 ( lim + ? ∞ → >> syms a b >> subs(a+b,a,4) ans = 4+b >> subs(cos(a)+sin(b),{a,b},{sym('alpha'),2}) ans = cos(alpha)+sin(2) >> syms x y >> f=x/y+y/x; >> [n,d]=numden(f) n = x^2+y^2 d = y*x 4 ? 符号微分: diff(S) 求符号表达式 S 对于默认自变量的微分 diff(S,v) 求符号表达式 S 对于自变量 v 的微分 diff(S,n) 求符号表达式 S 对于默认自变量的 n 次微分例:求axx44,的微分和三次微分 ? 符号积分 int(S) 求符号表达式 S 对于默认自变量的不定积分 int(S,v) 求符号表达式 S 对于自变量 v 的不定积分 int(S,a,b) 求符号表达式 S 对于默认自变量从 a 到b的定积分例: ∫ ∫ + + 1 0 2 ) 1 log( , ) 1 ( dx x x dz z x >> diff(x^4) ans = 4*x^3 >> diff(x^(4*a),x) ans = 4*x^(4*a)*a/x >> diff(x^4,2) ans = 12*x^2 >> syms x a >> limit(sin(x)/x) ans = 1 >> limit(1/x,x,0,'right') ans = inf >> v=(1+a/x)^x; >> limit(v,x,inf,'left') ans = exp(a) 5 ? 符号求和 symsum(S) 求符号表达式 S 对于默认自变量的不定和 symsum(S,v) 求符号表达式 S 对于自变量 v 的不定和 synsum(S,a,b) 求符号表达式 S 对于默认自变量从 a 到b的有限和 ? 泰勒级数展开 taylor(f) 计算符号表达式 f 在默认自变量等于 0 处的 5 阶Taylor 级数展开式 taylor(f,n,v) 计算符号表达式 f 在自变量 v=0 处的 n-1 阶Taylor 级数展开式 taylor(f,n,v,a) 计算符号表达式 f 在自变量 v=a 处的 n-1 阶Taylor 级数展开式 ? 多变元多项式系数求取 C = COEFFS(P) 返回多项式 P 的系数 C = COEFFS(P,X) 返回多项式 P 中变量 X 的系数 4) 符号方程求解 ? 符号代数方程求解 g=solve(eq) 求解符号表达式 eq=0 的代数方程,自变量为默认自变量 g=solve(eq,var) 求解符号表达式 eq=0 的代数方程,自变量为 var >> syms a b c x1 x2 x3 >> y=a*x1+b*x2; >> r=coeffs(y,x1) r = [ b*x2, a] >> r=coeffs(y,x3) r = a*x1+b*x2 >> [r,c]=coeffs(y,x1) r = [ b*x2, a] c = [ 1, x1] >> [r,c]=coeffs(y,x3) r = a*x1+b*x2 c = 1 >> y=a*x1; >> [r,c]=coeffs(y,x1) r = a c = x1 >> [r,c]=coeffs(y,x2) r = a*x1 c = 1 >> syms x z >> f=x/(1+z)^2; >> int(f) ans = 1/2*x^2/(1+z)^2 >> int(f,z) ans = -x/(1+z) >> f=x*log(1+x); >> int(f,0,1) ans = 1/4 6 g=solve(eq1,eq2,…,eqn,var1,var2,…,varn) 求解符号表达式 eq1,eq2,…,eqn 组成的代数方程组,自变量分别为 var1,var2,…,varn.方程组的解缺省存在结构变量中. 例1:求解代数方程 0 2 = + + c bx ax 和代数方程组 ? ? ? ? ? ? ? = + ? = ? + = + ? 0 4 2 0 5 10 2 2 z y x z y x z y x >> syms x y z a b c >> s1=a*x^2+b*x+c; >> s2=x^2-y^2+z-10; >> s3=x+y-5*z; >> s4=2*x-4*y+z; >> solve(s1) ans = [ 1/2/a*(-b+(b^2-4*a*c)^(1/2))] [ 1/2/a*(-b-(b^2-4*a*c)^(1/2))] >> [x y z]=solve(s2,s3,s4) x = [ -19/80+19/240*2409^(1/2)] [ -19/80-19/240*2409^(1/2)] y = [ -11/80+11/240*2409^(1/2)] [ -11/80-11/240*2409^(1/2)] z = [ -3/40+1/40*2409^(1/2)] [ -3/40-1/40*2409^(1/2)] >> double([x,y,z]) ans = 3.6481 2.1121 1.1520 -4.1231 -2.3871 -1.3020 >> eval([x,y,z]) ans = 3.6481 2.1121 1.1520 -4.1231 -2.3871 -1.3020 >> S=solve(s2,s3,s4); >> [S.x,S.y,S.z] ans = [ -19/80+19/240*2409^(1/2), -11/80+11/240*2409^(1/2), -3/40+1/40*2409^(1/2)] 7 [ -19/80-19/240*2409^(1/2), -11/80-11/240*2409^(1/2), -3/40-1/40*2409^(1/2)] FSOLVE 求解多元非线性方程组,其中方程的形式为:F(X)=0 X=FSOLVE(FUN,X0) 解初值为 X0 的由 FUN 表达的方程组. FUN 可以用@来表示 x = fsolve(@myfun,[2 3 4]) 其中 myfun 是如下 MATLAB 函数: function F = myfun(x) F = sin(x); 例: c = -1; % define parameter first x = fsolve(@(x) myfun(x,c),[-5;-5]) function F = myfun(x,c) F = [ 2*x(1) - x(2) - exp(c*x(1)) -x(1) + 2*x(2) - exp(c*x(2))]; ? 符号微分方程求解 r=dsolve('eq1,eq2,…','cond1,cond2,…','v')求由 eq1,eq2,…指定的常微分方程的符号解, 参数 cond1,cond2,…为指定常微分方程的边界条件或初始条件.自变量 v 如果不指定将为默认变量. 在方程中用大写字母 D 表示一次微分,D2 和D3 分别表示二次及三次微分,D 后面的字符为因变量. 例:求微分方程 ay dt dy = 的通解和当 y(0)=b 时的特解. 例:求微分方程 y a dt y d 2 2 2 ? = 当y(0)=1 及0π=??????adt dy 时的特解 > syms y a b >> dsolve('Dy=a*y') ans = C1*exp(a*t) >> dsolve('Dy=a*y','y(0)=b') ans = b*exp(a*t) 8 1.2. 数学模型之间的转换假设对象的传递函数为: n n n n n n n a s a s a s c s c s c s U s Y + + + + + + + = ? ? ? ? 1 1 1 2 2 1 1 ) ( ) ( ? ? 1.2.1 传递函数转换为可控标准型可控标准型为: u x x x x a a a a x x x x n n n n n n n ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 1 2 1 1 2 1 1 2 1 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? [ ] ? ? ? ? ? ? ? ? ? ? ? ? = = ? n n n x x x c c c CX y ? ? 2 1 1 1 状态变量的初值 X(0)求解方法: 0 0 0 0 0 1 2 1 1 2 1 x c x c x c x c y n n n n + + + + = ? ? ? ( ) n n n n n x x f x c x c x c x c y , , 1 1 2 3 1 2 1 ? ? ? ? = + + + + = ? ( ) ( ) ( ) 0 , , 0 ) 0 ( 1 1 n x x f y ? ? = ( ) 0 , , 0 0 1 2 n x x f y ? ? ? = 依次类推,可以得到一组线性代数方程组 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ? ? ? ? ? ? ? = = = ? ? 0 , 0 0 , , 0 ) 0 ( 0 , , 0 0 1 1 1 0 1 1 1 0 n n n n n x x f y x x f y x x f y ? ? ? ? ? 其中 ( ) ( ) ( ) n n c c c y y , , , ; , , 0 2 1 1 0 ? ? ? 均为已知,解方程可得到 0 , , 0 , 0 2 1 n x x x ? . >> syms y a >> dsolve('D2y=-a^2*y','y(0)=1,Dy(pi/a)=0') ans = cos(a*t) 9 假设系统的传递函数为 ( ) n n n n m m m m a s a s a s c s c s c s c s G + + + + + + + + = ? ? ? ? 1 1 1 1 1 1 0 ? ? ,则可控标准型为: u x x x x a a a a x x x x n n n n n n n ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 1 2 1 1 2 1 1 2 1 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? [ ] ? ? ? ? ? ? ? ? ? ? = ? n m m x x c c c c y ? ? ? 1 0 1 1 0 0 状态变量初值: n n x n q x q x q x n C x C x C y ) , 1 ( ... ) 2 , 1 ( ) 1 , 1 ( ) ( ... ) 2 ( ) 1 ( 2 1 2 1 + + + = + + + = [ ] u n q x n q x q x q u x n n A x n A x n A n q x n q x q x q x n q x q x q y n n n n ) , 1 ( ) , 2 ( ... ) 2 , 2 ( ) 1 , 2 ( ) , ( ... ) 2 , ( ) 1 , ( ) , 1 ( ) 1 , 1 ( ... ) 2 , 1 ( ) 1 , 1 ( ) , 1 ( ... ) 2 , 1 ( ) 1 , 1 ( 2 1 2 1 3 2 2 1 + + + + = + + + + + ? + + + = + + + = ? ? ? ? [ ] ? ? ? ? ? ? ? ? ? u n q u n q x n q x q x q u n q u x n n A x n A x n A n q x n q x q x q u n q x n q x q x q y n n n n ) , 2 ( ) , 1 ( ) , 3 ( ... ) 2 , 3 ( ) 1 , 3 ( ) , 1 ( ) , ( ... ) 2 , ( ) 1 , ( ) , 2 ( ) 1 , 2 ( ... ) 2 , 2 ( ) 1 , 2 ( ) , 1 ( ) , 2 ( ... ) 2 , 2 ( ) 1 , 2 ( 2 1 2 1 3 2 2 1 + + + + + = + + + + + + ? + + + = + + + + = u n j q u n q u n q x n j q x j q x j q y j j n j ) , ( ... ) , 2 ( ) , 1 ( ) , 1 ( ... ) 2 , 1 ( ) 1 , 1 ( ) 2 ( ) 1 ( 2 1 ) ( + + + + + + + + + + = ? ? 则可控标准型状态初值求取的方程组如下: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ) 0 ( ) 0 ( ) 0 ( 0 ) , 1 ( ) , 2 ( ) , 1 ( 0 0 0 ) , 1 ( ) , 2 ( 0 0 ) , 1 ( 0 0 0 0 0 ) 0 ( ) 0 ( ) 0 ( ) , ( ) 2 , ( ) 1 , ( ) , 2 ( ) 2 , 2 ( ) 1 , 2 ( ) , 1 ( ) 2 , 1 ( ) 1 , 1 ( ) 0 ( ) 0 ( ) 0 ( ) 2 ( 2 1 ) 1 ( n n n u u u n q n n q n n q n q n q n q x x x n n q n q n q n q q q n q q q y y y ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 例:求传递函数为 2 5 3 3 2 ) ( 2 3 2 + + + + + = s s s s s s G 的可控标准型实现,并求取状态初值. 10 解: 可控标准型实现: clear num=[1 2 3]; den=[1 3 5 2]; m=length(num); n=length(den)-1; A=[zeros(n-1,1),eye(n-1)]; A(n,:)=-den(n+1:-1:2)/den(1) B=zeros(n,1);B(n)=1 C=zeros(1,n);C(1:m)=num(m:-1:1)/den(1) 结果: A = 0 1 0 0 0 1 -2 -5 -3 B = 0 0 1 C = 3 2 1 状态初值求取(方法 1) y=[1 2 3];%y=[y(0) y'(0) y''(0)...] u=[0 4 5];%u=[0 u(0), u'(0),...] syms x1 x2 x3 xt xt=A(3,1)*x1+A(3,2)*x2+A(3,3)*x3+u(2); eq1=C(1)*x1+C(2)*x2+C(3)*x3-y(1); eq2=C(1)*x2+C(2)*x3+C(3)*(A(3,1)*x1+A(3,2)*x2+A(3,3)*x3+u(2))-y(2); eq3=C(3)*A(3,1)*x2+(C(1)+C(3)*A(3,2))*x3+(C(2)+C(3)*A(3,3))*(A(3,1)*x1+A(3,2)*x2+A(3,3)*x3+u(2))+C(3)*u(3 [x1,x2,x3]=solve(eq1,eq2,eq3); x0_1=double([x1 x2,x3]) 结果: x0_1 = -1 0 4 该方法完全最大的缺点是只能求解状态个数为 3 的系统,不通用. 状态初值求取(方法 2) syms x xt2 for i=1:n %确定状态变量 x(i)=['x' int2str(i)]; end eq(1)=x*C'; xt2=x*A(n,:)'; u_coef=zeros(n,n); eq_coef=zeros(n,n); for i=1:(n-1) %列写方程(只求取了前 n-1 个状态方程的系数) for j=1:n [ll,temp]=coeffs(eq(i),x(j)); if length(ll)==2 eq_coef(i,j)=ll(2); 11 elseif temp==x(j) eq_coef(i,j)=ll; end end eq(i+1)=x(2:n)*eq_coef(i,1:n-1)'+eq_coef(i,n)*xt2; end u_coef=zeros(n,n); for i=2:n %确定方程中由 u 和y各阶导数项组成的常数项 for j=2:i u_coef(i,j)=eq_coef(i-j+1,n); end end eq=eq-y+(u_coef*u')'; X0=solve(eq(1),eq(2),eq(3)); %解方程 x0_2=double([X0.x1 X0.x2 X0.x3]); 由以上程序可以看出,解方程时还不能做到通用,如何解决这个问题呢?我们采用 fsolve 函数, 状态初值求取(方法 3) 将方法二中的最后三行改写为: for j=1:n %求取最后一个状态方程的系数 ll=coeffs(eq(i),x(j)); if length(ll)==2 eq_coef(n,j)=ll(2); end end yu=-y+(u_coef*u')'; %常数项的系数 x0_3=fsolve(@(x) myfun(x,eq_coef,yu),[0 0 0]) function F = myfun(x,coef_x,coef) n=length(coef); for i=1:n F(i)=coef(i); for j=1:n F(i)=F(i)+x(j)*coef_x(i,j); end end 采用 fsolve 函数,简化上述编程过程,可以得到以下程序: 状态初值求取(方法 4) eq_coeff1=zeros(n,n); eq_coeff1(1,:)=C; for i=1:n-1 12 eq_coeff1(i+1,:)=[0 eq_coeff1(i,1:n-1)]+eq_coeff1(i,n)*A(n,:); end u_coef=zeros(n,n); for i=2:n for j=2:i u_coef(i,j)=eq_coef(i-j+1,n); end end yu=-y+(u_coef*u')'; x0_4=fsolve(@(x) myfun(x,eq_coef,yu),zeros(1,n)) function F = myfun(x,coef_x,coef) n=length(coef); for i=1:n F(i)=coef(i); for j=1:n F(i)=F(i)+x(j)*coef_x(i,j); end end 1.2.2 传递函数转换为可观标准型可观标准型为 Bu AX X + = ? CX y = ,其中: [ ] 1 0 0 , 1 0 0 0 1 0 0 0 1 0 0 0 1 2 1 1 2 1 ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? C c c c c B a a a a A n n n n n 由输入和输出的各阶导数的初值来计算状态变量初值 X(0)的矩阵表达式 13 ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( )( ) ( ) ( ) ( ) ( ) ( )? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 2 3 1 3 3 2 1 3 2 1 1 2 1 3 3 2 1 2 1 1 3 2 1 0 0 0 0 0 0 0 0 0 ) 0 ( ) 0 ( 0 0 0 1 0 1 1 0 0 0 0 0 n n n n n n n n n n n n n n n n n u u u u u c c c c c c c c y y y y y a a a a a a a x x x x x ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1.2.3 输入输出同阶次的处理方法假设对象传递函数为 n n n n n n n n a s a s a s c s c s c s c s U s Y s G + + + + + + + + = = ? ? ? ? 1 1 1 1 1 1 0 ) ( ) ( ) ( ? ? 则状态方程和输出方程可以表示为: u β β β β x x x x a a a a x x x x n n n n n n n n n ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1 2 1 1 2 1 1 2 1 1 2 1 1 0 0 0 1 0 0 0 0 1 0 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? [ ] u x x x y n 0 2 1 0 0 1 β + ? ? ? ? ? ? ? ? ? ? ? ? = ? ? 其中: ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? = ? = = ? ? 0 1 1 1 1 0 2 1 1 2 2 0 1 1 1 0 0 β a β a β a c β β a β a c β β a c β c β n n n n n ? ? 状态的初值可以通过下列方程获得: 14 ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? = ? ? = ? = ? ? ? ? ? ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( 1 2 ) 2 ( 1 ) 1 ( 0 ) 1 ( 2 1 0 3 1 0 2 0 1 u u u u y x u u u y x u u y x u y x n n n n n n β β β β β β β β β β ? ? ? ? ? ? ? ? ? ? 1.2.4 状态初值的通用求解方法状态空间模型一般由下式表示: DU CX Y BU AX X + = + = ? 根据状态与输入输出之间的关系,我们可以得到: ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 1 ( ) 2 ( ) 3 ( 2 1 ) 1 ( ) 1 ( ) 1 ( 2 ? ? ? ? ? ? ? ? + + + + + = + = + + + = + = + + = + = + = n n n n n n n n Du CBu CABu Bu CA X CA Du CX y u D u CB CABu X CA u D X C y u D CBu CAX u D X C y Du CX y ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 用矩阵表示: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( 0 0 0 ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 1 ( ) 3 ( 4 3 2 2 4 3 2 1 1 3 2 ) 1 ( ) 3 ( n n n n n n n u u u u u D CB B CA B CA B CA D CB CAB B CA D CB CAB D CB D x x x x x CA CA CA CA C y y y y y ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? )) 0 ( ) 0 ( ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( 1 U N Y Mo X U N X Mo Y * ? * = * + * = ? 式中,Mo 就是能观性判别矩阵,N 是一个下三角矩阵.通过此式就可以计算出状态初值. 实验内容 1.编写传递函数的可观标准型实现及其初值求取的通用程序,并以系统传递函数 ) 2 )( 3 )( 4 ( ) 1 ( 24 ) ( ) ( + + + + = s s s s s U s Y , 初始条件0)0(,5.0)0(,1)0(===yyy???,,1)0(=u15 0 ) 0 ( , 5 . 0 ) 0 ( = ? = u u ? ? ? 为例进行验证. 2.编写由传递函数转换为如下状态方程表达式以及状态初值求取的通用程序 u β β β β x x x x a a a a x x x x n n n n n n n n n ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1 2 1 1 2 1 1 2 1 1 2 1 1 0 0 0 1 0 0 0 0 1 0 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? [ ] ? ? ? ? ? ? ? ? ? ? ? ? = n x x x y ? ? 2 1 0 0 1 并以系统传递函数 ) 2 )( 3 )( 4 ( ) 1 ( 24 ) ( ) ( + + + + = s s s s s U s Y ,初始条件 0 ) 0 ( , 5 . 0 ) 0 ( , 1 ) 0 ( = = = y y y ? ? ? , 0 ) 0 ( , 5 . 0 ) 0 ( , 1 ) 0 ( = ? = = u u u ? ? ? 为例进行验证. 16 实验二连续系统离散化实验目的和要求: 1) 熟悉Matlab矩阵引用及高维矩阵操作 2) 实现4 阶Runge_Kutta常微分方程(组)数值解 3) 应用Runge_Kutta方法解常微分方程(组) 实验方法: 2.1 Matlab 基本函数 (1) 字符串运算: EVAL 执行 MATLAB 表示的字符串. EVAL(s), 执行字符串 s 表示的运算. 例: >> t=1; >> x='t^2+t+1'; >> c=eval(x) c = 3 >> syms a s d >> a=s+d; >> s=1; >> d=2; >> f=eval(a) f = 3 (2) ODE45 求解非病态或非刚性(non-stiff)微分方程使用句法: [T,Y] = ode45(odefun,tspan,y0) [T,Y] = ode45 (odefun,tspan,y0,options) [T,Y,TE,YE,IE] = ode45 (odefun,tspan,y0,options) sol = solver(odefun,[t0 tf],y0...) 输入输出项的描叙: Odefun:需要积分的函数,先把要积分的函数 ( , ) y f t y = ? 写成 XXX.m 文件的格式.调用时函数名前面加上符号:@; Tspan: 积分区间.写成 0 [ , ] f t t 的形式, 0 t -积分起始时间, f t -结束时间; Y0: 设定一个积分初值,默认为 0; Options:积分选项,描叙积分的属性,如积分精度等,一般取默认值.也可以用 odeset 进行设置. T: 积分点.是一个向量. Y: 对应于积分点处的积分值. 17 例: 求解如下方程: 2 1 3 3 1 2 3 2 1 51 . 0 y y y y y y y y y ? = ? = = ? ? ? 1 ) 0 ( 1 ) 0 ( 0 ) 0 ( 3 2 1 = = = y y y 步骤如下: 1、建立积分函数的 M 文件: function dy = rigid(t,y) dy = zeros(3,1); % a column vector dy(1) = y(2) * y(3); dy(2) = -y(1) * y(3); dy(3) = -0.51 * y(1) * y(2); 2、再在命令窗口(command window)或M文件中写如下: options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]); [T,Y] = ode45(@rigid,[0 12],[0 1 1],options); 3、做图 plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.') 得到图形如下: 注:ode45 函数只是作为与自己编程结果的比较,本次实验不能采用 ode45 来完成实验作业内容. 2.2 欧拉法假设原始微分方程为 ) , ( t y f y = ? ,欧拉法有四种形式: ? 前差公式: n n n hf y y + = +1 18 ? 后差公式: 1 1 + + + = n n n hf y y ? 梯形积分公式(平均公式) : ) ( 2 1 1 + + + + = n n n n f f h y y ? 中心差公式: n n n hf y y 2 1 1 + = ? + 欧拉前差方法的 MATLAB 实现例:用欧拉前差方法求解一阶微分方程 1 . 0 , 1 ) 0 ( ] 1 , 0 [ 2 = = = ? = h y t y y ? function [yt,th]=oulaqiancha(f,y0,t0,td,h) %format long; y=y0; t=t0; yt(1)=double(y0); th(1)=double(t0); for i=2:fix((td-t0)/h)+1 fn=double(eval(f)); yt(i)=yt(i-1)+fn*h; th(i)=th(i-1)+h; y=yt(i);t=th(i); end >> [yt,t]=oulaqiancha('-y^2',1,0,1,0.1) yt = 1.0000 0.9000 0.8190 0.7519 0.6954 0.6470 0.6052 0.5685 0.5362 0.5075 0.4817 t = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 例:用欧拉前差方法求解一阶微分方程组 ) , , ( ) , , ( 2 1 2 2 2 1 1 1 t y y f y t y y f y = = ? ? 解:方法一,将f和y表示成向量的形式 ) ), 2 ( ), 1 ( ( ) 2 ( ) ), 2 ( ), 1 ( ( ) 1 ( 2 1 t y y f y t y y f y = = ? ? 则欧拉前差法的子程序为: function [yt,th]=oulaqiancha_1(f,y0,t0,td,h) %format long; y=y0; t=t0; yt(:,1)=double(y0); th(1)=double(t0); 19 for i=2:fix((td-t0)/h)+1 fn=double(eval(f)); yt(:,i)=yt(:,i-1)+fn'.*h; th(i)=th(i-1)+h; y=yt(:,i);t=th(i); end 主程序 >>syms f >> f(1)='-y(1)^2'; >> f(2)='-y(2)^2'; >> [yt,t]=oulaqiancha_1(f,[1 1]',0,1,0.1); >> f(1)='-y(1)^2+y(2)'; >> f(2)='-y(2)^2'; >> [yt,t]=oulaqiancha_1(f,[1 1]',0,1,0.1); 注:该子程序也可以求解单个微分方程. 方法二:将f用matlab 中函数的形式表示 function dy=oulaqiancha_fun(y,t) dy(1)=-y(1)^2+y(2); dy(2)=-y(2)^2+t; dy=dy'; function [yt,th]=oulaqiancha_2(f,y0,t0,td,h) %format long; yt(:,1)=double(y0); th(1)=double(t0); for i=2:fix((td-t0)/h)+1 fn=feval(f,yt(:,i-1),th(i-1)); yt(:,i)=yt(:,i-1)+fn*h; th(i)=th(i-1)+h; y=yt(:,i);t=th(i); end 主程序: >> clear >> [yt,t]=oulaqiancha_2('oulaqiancha_fun',[1 1]',0,1,0.1) yt = 1.0000 1.0000 0.9900 0.9749 0.9579 0.9411 0.9258 0.9131 0.9033 0.8969 0.8940 1.0000 0.9000 0.8290 0.7803 0.7494 0.7332 0.7295 0.7363 0.7521 0.7755 0.8054 t = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 20 2.3 龙格-库塔法 (1) 二阶龙格-库塔法: 二阶龙格-库塔公式,也称为改进欧拉公式. ) 2 1 , 2 1 ( 1 n n n n n hf y h t hf y y + + + = + ? ? ? ? ? ? ? ? ? ? ? ? + + + + = + n n n n n n hf y h t f f h y y 3 2 , 3 2 3 4 1 [ ] ) , ( 2 1 n n n n n n hf y h t f f h y y + + + + = + (2) 三阶龙格-库塔公式: ) 3 2 , 3 2 ( ) 3 , 3 ( ) , ( ] 3 [ 4 2 3 1 2 1 3 1 1 k h y h t f k k h y h t f k y t f k k k h y y n n n n n n n n + + = + + = = + + = + (3) 四阶龙格-库塔公式: ) , ( ) 2 , 2 ( ) 2 , 2 ( ) , ( ] 2 2 [ 6 3 4 2 3 1 2 1 4 3 2 1 1 hk y h t f k k h y h t f k k h y h t f k y t f k k k k k h y y n n n n n n n n n n + + = + + = + + = = + + + + + = + 四阶龙格-库塔算法实现 function [T]=Runge_Kutta(f,x0,y0,h,n) % [T]=Runge_Kutta(f,x0,y0,h,n) % f 待解方程(组) % x0 初试自变量值 % y0 初试函数值 % h 步长 % n 步数 if nargin<5 n=100; end r=size(y0); r=r(1); s=size(x0); s=s(1); 21 r=r+s; T=zeros(r,n+1); T(:,1)=[y0;x0]; for t=2:n+1 k1=feval(f,T(1:r-1,t-1)); k2=feval(f,[k1*(h/2)+T(1:r-1,t-1);x0+h/2]); k3=feval(f,[k2*(h/2)+T(1:r-1,t-1);x0+h/2]); k4=feval(f,[k3*h+T(1:r-1,t-1);x0+h]); x0=x0+h; T(:,t)=[T(1:r-1,t-1)+(k1+k2*2+k3*2+k4)*(h/6);x0]; end 例:求解下列著名的洛伦兹吸引子微分方程组 3 2 1 3 2 1 3 2 1 2 1 ) ( ) ( Bx x x dx x x x R dx x x Q dx ? = ? ? = ? = 其中Q,B,R 为常数. 步骤: 1) 新建函数fun.m function dy=fun(x) R=28.025; B=8.0/3; Q=10.0; dy(1)=Q*(x(2)-x(1)); dy(2)=((R-x(3))*x(1)-x(2)); dy(3)=(x(1)*x(2)-B*x(3)); dy=dy'; 2) 运算: >> T=Runge_Kutta('fun',0,[10.1;9.8;11.2],0.01,5000); >> plot3(T(1,:),T(2,:),T(3,:)) 3)观察图形窗口, 改变角度, 如图 22 2.4 亚当斯-巴什福思(Adams-Bashforth) (1)亚当斯-巴什福思显式公式(Adams-Bashforth) 由泰勒级数展开式推导出来,其中 n y ? ? 由后项差分代替. [ ] ) ( 3 2 3 1 1 h O f f h y y n n n n + ? + = ? + 求解一个新的 y 值,需要 f 的两个值.由于在零点处只有 y0 和f,因此不能自启动,一般用同阶(二阶)龙格库塔法来启动. (2)亚当斯-莫尔顿隐式公式(Adams-moulton) 由向后展开的泰勒级数公式推导得到的,截掉 fn+1 以后各项得到 ) ( 2 1 1 h O hf y y n n n + + = + + 隐式:因为 yn+1 中包含有 fn+1,而fn+1 反过来又包含 yn+1. 步骤:先估算一个 yn+1,计算 fn+1,然后用隐式公式求得 yn+1 的新估计值,重复迭代, 直至前后两次的 yn+1 值之间的差在要求的范围内. 2.5 预估-校正法使用隐式公式必须包含一个能为每步解都提供一个首次精确估值的方法,一般选择误差阶数与隐式相同的显式公式. (1)选择四阶 Adams 显示公式作为"预估" ,四阶 Adams 隐式公式作为"校正" ,计算过程采用四阶龙格-库塔公式开始,以便得到初始条件之外的最初三步 h 的y和f值. 四阶 Adams 显式公式: ? ? ? ? ? ? ? + ? + = ? ? ? + 3 2 1 0 1 24 9 24 37 24 59 24 55 n n n n n n f f f f h y y 四阶 Adams 隐式公式: ? ? ? ? ? ? + ? + + = ? ? + + + 2 1 1 1 1 24 1 24 5 24 19 24 9 n n n i n n i n f f f f h y y (2)汉明(Hamming)积分法 ? 预估公式: ( ) 2 1 3 ) 0 ( 1 2 2 3 4 ? ? ? + + ? + = n n n n n f f f h y y ,不考虑稳定性 ? 修正公式: ) ( 121 112 ~ ) 0 ( ) 0 ( 1 ) 0 ( 1 n n n n y y y y ? + = + + ,改善预估值,降低校正迭代次数.稳定性要好. 23 ? 校正公式: ) 2 ( 8 3 ) 9 ( 8 1 1 ) ( 1 2 ) 1 ( 1 ? + ? + + ? + + ? = n n i n n n i n f f f h y y y 其中: ) , ( )) ( , ( n n n n n y t f t y t f f = = ) , ( )) ( , ( 1 1 1 1 1 ? ? ? ? ? = = n n n n n y t f t y t f f ) , ( )) ( , ( 2 2 2 2 2 ? ? ? ? ? = = n n n n n y t f t y t f f ) , ( )) ( , ( 3 3 3 3 3 ? ? ? ? ? = = n n n n n y t f t y t f f ) , ( ) , ( ) ( 1 1 ) ( 1 ) ( ) 1 ( 1 i n n i t n i n y t f y t f f n + + + + = = + 2.6 状态转移法离散化假设对象的状态方程: BU AX X + = ? ,则离散化之后的状态方程为: [ ] ) ( ) ( ) ( ) ( ) 1 ( kT U T kT X T T k X m Φ + Φ = + AT e T = Φ ) ( ∫ ∫ ? = ? Φ = Φ T T A T m Bd e Bd T T 0 ) ( 0 ) ( ) ( τ τ τ τ 矩阵指数的 matlab 函数: EXPM 矩阵指数例:若???????=1100A,则 ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? + + = ? ? ? ? ? ? + ? = ? = = Φ ? ? ? ? ? ? T T AT e e s s s s L s s L A sI L e T 1 0 1 1 1 ) 1 ( 1 0 1 1 1 0 ) ( ) ( 1 1 1 1 matlab 编程如下: >> syms T >> A=[0 0;1 -1]; >> expm(A*T) ans = [ 1, 0] [ -exp(-T)+1, exp(-T)] 结果一致. 2.7 连续系统传递函数到离散系统模型的转化 C2D 将连续时间模型转化为离散时间模型 SYSD = C2D(SYSC,Ts,METHOD)——将连续模型转换为离散模型, METHOD 缺省为采 24 用零阶保持器的方法,Ts 为采样周期. Method: 'zoh'——采用零阶保持器 'foh'——采用一阶保持器 'tustin'——采用双线形(tustin)逼近方法 'prewarp'——采用改进的 tustin 方法 'matched'——采用 SISO 系统的零极点匹配法. 对于状态方程来说, [SYSD,G] = C2D(SYSC,Ts,METHOD) 返回的G 表示连续时间系统的初始条件与离散时间系统初始条件之间的转换系数. 假设 x0,u0 是SYSC 的初始状态和初始输入, 则SYSD 中等价的初始条件为 xd[0] = G * [x0;u0], ud[0] = u0 [numd,dend]=tfdata(sysd,'v') 返回传递函数 sysd 中分子和分母的系数除了采用以上 Matlab 中函数 c2d 进行离散化之外,还可以根据自动控制原理中介绍的 Z 变换方法来对连续系统进行离散化.首先复习拉氏变换和 Z 变换. Laplace 变换和 Ztrans 变换 L = LAPLACE(F) 系统 F(缺省独立变量为 t)的拉氏变换,缺省返回为 s 的函数,如果 F = F(s),那么 LAPLACE 返回为 t 的函数 L = L(t). F = ILAPLACE(L) 系统 L(缺省独立变量为 s)的拉氏逆变换,缺省返回为 t 的函数,如果 L = L(t),那么 LAPLACE 返回为 s 的函数 L = L(s). F = ZTRANS(f) 表示函数 f(n)的Z变换,返回为 z 的函数 F(z) ,f = f(n) => F = F(z). F(z) = symsum(f(n)/z^n, n, 0, inf) (离散时间函数) f = IZTRANS(F) ——表示 F(z)的反变换,缺省时返回 f 为n的函数 f(n),F = F(z) => f = f(n). 如果 F = F(n), 那么返回 f 为k的函数: f = f(k). 假设系统的传递函数为 G(s),则采用零阶保持器的 Z 变换为: ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? s s G Z z s G s e Z Ts ) ( ) 1 ( ) ( 1 1 由于 laplace 变换和 ztrans 变换都是从时域出发的,因此,在已知对象传递函数的前提下,首先需要采用 laplace 反变换将其转化为连续时域形式,然后进行采样(采样周期为 T)转化为离散时域形式,最后再采用 Z 变换的方法将其转换了离散系统模型.具体流程如下图所示: 25 Matlab 编程如下: clear num=1; den=[1 1 0]; syms w n T z wnum=poly2sym(num, 'w '); %求以符号 w 表示的传递函数 G(w) wden=poly2sym(den, 'w'); sys=wnum/wden/w; %G(w)/w ilsys=ilaplace(sys); %拉氏反变换 ilsys_T=subs(ilsys,'t',n*T); %离散时间表示 zsys_1=ztrans(ilsys_T); zsys_2=(z-1)/z; zsys=zsys_1*zsys_2; sysd=simplify(zsys); [numd dend]=numden(sysd); Ts=1; numd1=sym2poly(subs(numd,'T',Ts)); dend1=sym2poly(subs(dend,'T',Ts)); ml=max(length(numd1),length(dend1)); numd1=[zeros(1,ml-length(numd1)) numd1] numd2=[zeros(1,ml-length(dend1)) dend1] 运行结果: numd1 = 0 0.3679 0.2642 numd2 = 1.0000 -1.3679 0.3679 直接采用 c2d 的方法 >> clear >> num=1; >> den=[1 1 0]; >> sysc=tf(num,den); >> sysd=c2d(sysc,1) Transfer function: 传递函数 G(s)/s 拉氏反变换连续脉冲响应模型g1(t) 采样 t——〉nT 离散脉冲响应模型g1(nT) Z变换传递函数G(z) 差分方程系数考虑零阶保持 (z-1)/z 考虑零阶保持 26 0.3679 z + 0.2642 z^2 - 1.368 z + 0.3679 Sampling time: 1 实验内容 1. 编写欧拉后差公式, 计算实验 1 中习题 2 所得状态方程的结果 (其中输入 u=1, 步长 h=0.02, ] 1 , 0 [ ∈ t ) ,并绘制状态和输出曲线. 2.编写二阶、三阶龙格库塔法通用程序,然后以实验 1 中习题 2 所得状态方程为例,分别采用二阶、三阶龙格库塔法进行求解(其中输入 u=1,步长 h=0.02, ] 1 , 0 [ ∈ t ) ,并与四阶龙格库塔法结果进行比较.绘制状态与输出曲线. 3.编写汉明积分法通用程序,并用例 3-5 进行验证. 4.编写用状态转移法对连续系统状态方程进行离散化的通用程序.并计算实验 1 中习题 2 所得状态方程进行离散化. 27 实验三连续控制系统仿真实验目的和要求: 掌握控制系统仿真的两种方法:与信息流方向无关的仿真方法和与信息流方向有关的仿真方法; 掌握非线性环节的 Matlab 实现; 掌握 PID 控制系统仿真方法. 实验方法: 3.1 面向闭环系统传递函数的数字仿真 ? 先写出闭环系统传递函数 ? 用Z变换方法获得闭环系统的差分方程迭代格式或者称仿真模型; ? 在计算机上进行仿真 3.2 面向系统结构的数字仿真常用的方法: 按环节分别建立离散化模型, 然后引入反馈通道的表达式来构成仿真模型. 3.2.1 与信息流无关的仿真方法仿真计算过程为:按环节的编号顺序依次计算,每计算一个动态响应点,各环节都要计算一次,待各环节计算完后,再对各环节的输入信号进行综合,以便得到下一次计算时各环节差分方程的输入量.以下列控制系统为例进行介绍环节一: 1 ) ( ) ( ) ( 1 2 1 ? = = z T z E z E z G 差分方程: ) 1 ( ) 1 ( ) ( 1 2 2 ? + ? = n TE n E n E 环节二: T T e z e z E z Y z G ? ? ? ? = = 1 ) ( ) ( ) ( 2 2 差分方程: ) 1 ( ) 1 ( ) 1 ( ) ( 2 ? ? + ? = ? ? n E e n y e n y T T 比较环节: ) 1 ( ) ( ) ( 1 ? ? = n y n R n E matlab 编程如下: clear; r=1;tfinal=200;T=0.1; G1(z) R(z) Y(z) G2(z) T z-1 E1(z) T E2(z) 28 sp=r*ones(tfinal,1); e1_0=0;e2_0=0;y_0=0; %仿真计算顺序与信息流无关 for n=1:tfinal e1(n)=sp(n)-y_0; e2(n)=e2_0+T*e1_0; y(n)=exp(-T)*y_0+(1-exp(-T))*e2_0; e1_0=e1(n); e2_0=e2(n); y_0=y(n); end figure(1); plot([1:tfinal],y,'k',[1:tfinal],sp,'k'); 0 50 100 150 200 0 0.2 0.4 0.6 0.8 1 1.2 1.4 对于复杂的控制系统,可以将控制系统的结构输入到计算机中,控制系统的结构可以用连接矩阵来表示各个环节之间的关系. ? ? ? = + = CX Y PX FM U 其中:U——环节输入, M——系统输入, X——环节输出, Y——系统输出 F——它表示 r 个外部输入量与各环节输入端的连接关系,称为输入连接矩阵.是n*r 矩阵; P——它表示各环节输出量与系统各环节输入端的连接关系,称为环节连接矩阵.是n*n 矩阵; C——它表示系统 l 个输出量与各环节输出端的连接关系,称为输出连接矩阵.是l*n 矩阵; 仿真步骤如下图所示 29 3.2.2 与信息流相关的仿真方法要求在仿真过程中计算顺序按照信息流传递的环节次序进行.对于上一小节中的例子, 采用与信息流相关的仿真方法,则其对应的差分方程为: 环节一: 1 ) ( ) ( ) ( 1 2 1 ? = = z T z E z E z G 差分方程: ) ( ) 1 ( ) ( 1 2 2 n TE n E n E + ? = 环节二: T T e z e z E z Y z G ? ? ? ? = = 1 ) ( ) ( ) ( 2 2 差分方程: ) ( ) 1 ( ) 1 ( ) ( 2 n E e n y e n y T T ? ? ? + ? = 比较环节: ) 1 ( ) ( ) ( 1 ? ? = n y n R n E matlab 编程如下: clear; r=1;tfinal=200;T=0.1; sp=r*ones(tfinal,1); %仿真计算顺序与信息流有关 e11_0=0;e21_0=0;y1_0=0; for n=1:tfinal e11(n)=sp(n)-y1_0; e21(n)=e21_0+T*e11(n); y1(n)=exp(-T)*y1_0+(1-exp(-T))*e21(n); e21_0=e21(n); 参数初始化 (连接矩阵,参数初值) 各环节离散系数计算用连接矩阵综合信号,环节输入ui 计算环节输出(ui?xi) 计算系统输出y=CX 30 y1_0=y1(n); end figure(2); plot([1:tfinal],y1,'k',[1:tfinal],y,'k:',[1:tfinal],sp,'k'); 0 50 100 150 200 0 0.2 0.4 0.6 0.8 1 1.2 1.4 3.3 非线性环节仿真非线性环节: 饱和非线性、分段线性非线性、失灵区非线性、间歇非线性、继电非线性、具有死区的继电非线性、具有滞环的继电非线性、具有死区和滞环的继电非线性. ? 饱和、分段线性及失灵区非线性特性 u 出(m) u 入(m) h0 h1 -c1 c1 0 当1)(cmu<入时, ) ( ) ( 0 m u h m u 入出=当1)(cmu≥入时,若1)(cmu≥入,则 [ ] 1 1 1 0 ) ( ) ( c m u h c h m u ? + = 入出,否则 [ ] { } 1 1 1 0 ) ( ) ( c m u h c h m u ? ? + ? = 入出h1=0,为饱和非线性 h0=0,失灵区非线性 ? 间歇非线性特性 31 u 出(m) u 入(m) c1 -c1 45° 0 a b c d 在bd 段,输入 u 增加时,当前输出值 u 出(m)总大于前一时刻的输出值 u0 出(m) 在ca 段,输入 u 减小时,当前输出值 u 出(m)总小于前一时刻的输出值 u0 出(m) 在其他段,输出保持不变. 所以:设)(0mu入为前一次输入, ) ( 0 m u出为前一次输出,则间隙非线性可表示为当,)()(), ( ) ( 1 0 0 c m u m u m u m u ? < > 入出入入且则,)()(1cmumu?=入出即系统工作在右边的特性上. 当100)(), ( ) ( c m u u m u m u + > < 入出入入且,则 1 ) ( ) ( c m u m u + = 入出即系统工作在左边的特性上. 对于其他情况 ) ( ) ( 0 m u m u 出出=?继电及具有死区的继电非线性特性. u 出(m) u 入(m) c1 -c1 h1 -h1 0 u 出(m) u 入(m) c1 -c1 具有死区的继电非线性,当死区为 0 时,即为继电非线性. ? 具有滞环的继电非线性 32 c1 -c1 h1 -h1 u 出(m) u 入(m) 当输入递增且输入大于滞环区(h1)时,输出=c1 当输入递减且输入小于滞环区(-h1)时,输出=-c1 其余情况下,输出等于前一时刻的输出.即: 当110)(,)(), ( ) ( c m u h m u m u m u = ≥ > 出入入入则且;当110)(,)(), ( ) ( c m u h m u m u m u ? = ? ≤ < 出入入入则且;对于其他情况 ) ( ) ( 0 m u m u 出出=?具有死区和滞环的继电非线性特性 u 出(m) u 出(m) c1 -c1 -h1 h1 h0 -h0 0 当0)(hmu≤入时, 0 ) ( = m u出当1)(hmu≥入时, c m u = ) ( 出当1)(hmu?≤入时, c m u ? = ) ( 出其余情况,即0110)(,)(hmuhhmuh?<t1>τ; (iii)计算式(6-4)得K,按式(6-5)将)(1ty、)(2ty归一化,得)(1ty、)(2ty;(iv)按式(6-8)和式(6-9)计算 T 和τ; 上面 i~iv 已完成了一阶加纯滞后对象的建模.模型是否满足要求可通过仿真进行验证,并与实验曲线比较. 在响应曲线上根据时刻 t 读取响应数据或根据 y 的值获得对应的时刻可以用函数 interp1(),求解方程组可以用函数 fsolve(). INTERP1 为一维插值函数,对数据点之间函数的进行估值,这些数据点是由某些集合给定,具体用法如下: (i) yi=interp1(x,y,xi) 对数据点(x,y)进行插值,计算插值点 xi 的函数值.x 为维数为 N 的向量,y 为对应的函 38 数值.如果 y 为矩阵,则插值对 y 的每一列进行,若y的维数超出 x 或xi 的维数,则返回 NaN. (ii) yi=interp1(y,xi) 此格式默认 x=1:n,n 为向量 y 的元素个数值,或等于矩阵 y 的size(y,1). (iii) yi=interp1(x,y,xi,method) 此格式用 method 来指定插值的算法.默认为线性算法.其值常用的可以是如下的字符串: nearest——线性最近项插值. linear——线性插值. spline——三次样条插值. cubic——三次插值. (iv) yi=interp1(x,y,xi,mehod,'extrap') 利用指定的插值算法对不处于 x 区间内的 xi 进行外插. (v) pp=interp1(x,y,method,'pp') 利用指定的插值算法产生 y 的分段多项式 pp,pp 可以通过 ppval 估计点 xi 处的函数值. 所有的插值方法要求 x 是单调的,但x可以非连续等距. 例:根据比较粗糙的正弦曲线插值得到更精细的结果. x = 0:10; y = sin(x); %已知数据点 xi = 0:.25:10; %期望计算的点 yi = interp1(x,y,xi); %采用线性插值 plot(x,y,'o',xi,yi) %已知数据点用圆圈表示,插值得到的结果用线表示得到的结果如下: 0 1 2 3 4 5 6 7 8 9 10 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 在这个例子中, 函数 interp1 可用来在任何 x 的值上估计 y 值, 但不能找出对应于 y 值的 x 值. fsolve 的用法在实验一中有详细介绍,在此不赘述. (2) 建立二阶或 n 阶惯性环节传递函数模型的步骤如下: (i) 输入 x0,并根据阶跃响应曲线确定新稳态 y(∞); (ii) 根据式(6-4)和(6-5)计算 K 和y(t)的归一化值 y1(t); (iii) 在响应曲线上找出使 y1(t1)和y1(t2)分别为 0.4 和0.8 的t1 和t2; 39 (iv) 计算比值λ=t1/t2,如λ>0.46,则转(vii); (v) 如λ≤0.32,转(viii); (vi) 解联立方程(6-16),得到 T1 和T2,转(ix); (vii) 给出λ的信息,或是查表 6-1 得到 n,再由式(6-21)计算 T,求得如式(6-20)所示的传递函数,转(ix);或是退出计算,准备改用它法; (viii)给出λ提示信息,并按一阶传递函数拟合; (ix) 至此已完成了二阶加纯滞后对象或 n 阶惯性环节的建模.然而,所建立的模型是否满足要求则要通过仿真进行验证,并与实验曲线进行比较. 例4.3 书中习题 2,已知一温度对象的阶跃响应实验结果如下: t(s) 0 10 20 30 40 50 60 70 80 90 100 150 θ(℃) 0 0.16 0.65 1.15 1.52 1.75 1.88 1.94 1.97 1.99 2.00 2.00 阶跃扰动量为 Δq=1t/h.试采用二阶或 n 阶惯性环节建立该对象的传递函数. 解:此处给出一种方法供参考.首先根据具体的工业现场数据,画出拟合曲线,并且利用数据求得对象模型,画图表示,和习题 1 是一样的.主程序如下: %主程序 %现场数据如下 tin=[0:10:100 150]; ui=1; yout=[0 0.16 0.65 1.15 1.52 1.75 1.88 1.94 1.97 1.99 2.00 2.00]; GetModel (tin,ui,yout,0,'linear'); 其中调用的函数分别为: getorder.m 为根据数据,给出对象最适合的仿真模型的阶数: 40 function kn=getorder(kt) % 根据输入的比值kt(t1/t2)确定系统的阶次kn format short ff=[0.32 0.46 0.53 0.58 0.62 0.65 0.67 0.685 0.6975 0.71 0.7225 0.735 0.7425 0.75]; %阶次与比值关系 kn=0; if kt<=ff(1) kn=1; else if kt>=ff(length(ff)) kn=length(ff); end end for i=1:13 if 2*ktY(n) || ye if y>ytemp(j) min=i; i=floor((min+max)/2); j=j+1; 41 ytemp(j)=Y(i); else max=i; i=floor((min+max)/2); j=j+1; ytemp(j)=Y(i); end; if j>2 && (ytemp(j)==ytemp(j-1)||ytemp(j)==ytemp(j-2)) disp('step is too big, cannot find x in function GetXGivenY'); return; end end x=X(i); GetModel.m 为根据阶跃响应数据获得对象的模型. function GetModel(tin,ui,yout,n,method) % 根据具体的工业现场数据,求出对象模型,并以图形表示 % 参数 % tin:采样时间序列 % uin:阶跃输入幅值 % yout:采样时间点上的输出序列 % n:模型类型 % n=1——一阶纯滞后 % n=2——二阶惯性环节 % n=0——根据输入的数据选择合适的阶次 % method: 'linear':线性插值 % 'cubic':三次插值 % 'nearest':最近点插值 % 'spline':三次样条插值 % initializing... format long len=length(tin); uin=ui*ones(1,len); ti=tin(1):0.005:tin(len); %以0.005为步长插值中间点 yi=interp1(tin,yout,ti,method); e=1e-3; %二分法的误差 % plot the input and the step response; plot([tin(1) tin],[0 uin],'g',ti,yi,'r');grid;hold; % calculate y (∞) yf=max(yi)-yout(1); y1=yi/yf; %归一化 switch n case 1 % first-order model with lag 42 % calculate { t1,y(t1)},{t2,y(t2)} a1=1-exp(-0.5); a2=1-exp(-1); t1=GetXGivenY(a1, y1, ti, e); t2=GetXGivenY(a2, y1, ti, e); % calculate t he pa rameters: K=yf/ui T=2*(t2-t1) tao=2*t1-t2 % check t he pa rameters: sys=tf([K],[T 1],'inputdelay',tao) step(sys*ui,tin(len));grid; xaxis=tin(1)-tin(len)/10; yaxis=max(max(uin),max(yi))*1.35; axis([xaxis tin(len) min(yout) yaxis]); legend('Input','Step Response','Simulate Model Response'); title('Graph of 1st(lag) Model Response'); case 2 % second-order m odel % calculate { t1,y(t1)},{t2,y(t2)} a1=0.4;a2=0.8; t1=GetXGivenY(a1, y1, ti, e); t2=GetXGivenY(a2, y1, ti, e); % calculate t he pa rameters: K=yf/ui if t1/t2>0.46 disp(the model n is not fitalbe! end [T1 T2]=solve('T1+T2=(t1+t2)/2.16','T1*T2/(T1+T2)^2=1.74*t1/t2-0.55','T1,T2'); T1=eval(T1(1)) T2=eval(T2(1)) sys=tf([K],[(1.74*t1/t2-0.55)*((t1+t2)/2.16)^2 (t1+t2)/2.16 1]) step(sys*ui,tin(len));grid; xaxis=tin(1)-tin(len)/10; yaxis=max(max(uin),max(yi))*1.35; axis([xaxis tin(len) min(yout) yaxis]); legend('Input','Step Response','Simulate Model Response'); title('Graph of 2nd Model Response'); case 0 % the order n is selected by the system according to the data % calculate { t1,y(t1)},{t2,y(t2)} a1=0.4;a2=0.8; t1=GetXGivenY(a1, y1, ti, e); t2=GetXGivenY(a2, y1, ti, e); 43 % calculate t he pa rameters: K=yf/ui N=getorder(t1/t2) T=(t1+t2)/2.16/N sys=zpk([],[-1/T*ones(1,N)],K)/(T^N) % plot step(ui*sys,tin(len));grid; xaxis=tin(1)-tin(len)/10; yaxis=max(max(uin),max(yi))*1.35; axis([xaxis tin(len) min(yout) yaxis]); legend('Input','Step Response','Simulate Model Response'); title('Graph of Nth Model Response'); otherwise % wrong model type disp(The model type n is wrong! end 运行程序得到的结果如下: K = 2 N = 3 T = 10.1590 Zero/pole/gain: 0.0019076/(s+0.09844)^3 0 50 100 150 0 0.5 1 1.5 2 2.5 Graph of Nth Model Response Time (sec) Amplitude Input Step Response Simulate Model Response 图中的拟合三阶系统仿真后的曲线几乎与输入的数据得到的曲线完全重合. 4.1.3 从响应曲线求对象的微分方程模型 (1) 阶跃响应法如果已知自衡对象的阶跃响应曲线,利用书中式(6-25)~(6-30),可以得到对象微分方程模型的系数.其中阶乘函数为 factorial,求积分可以采用针对离散点数据积分的梯形积分函数trapz. 44 根据阶跃响应曲线求微分方程的例程如下,函数返回微分方程的系数: function a=GetDEFromStepResponse(x0, y, t, e) % pram\ x0: magnitude of step input % pram\ y: response of step input % pram\ t: time v ector % pram\ e: error t hreshold % get len, yEnd, a(1) len=length(y); yEnd=y(len); a(1)=x0/yEnd; tr=zeros(1, len); k=1; while true sum=0; for i=0:k-1 sum=sum+(-1)^(k+i+1)*a(i+1)/a(1)*(t.^(k-i-1))/factorial(k-i-1); end tr=sum.*(yEnd-y); aTemp=a(1)/yEnd*trapz(t,tr); aMax=max(a); if aMax~=0&&aTemp/aMax1 zz(1:N)=0; for i=1:N zz(i)=1-wi(i,L)/wi(i,L-1); end z=mean(abs(zz)); end if z<=sigma num=[fliplr(b) b0]; 50 den=[fliplr(a) 1]; sys=tf(num,den) figure; plot(P,Q,'rx'); hold on nyquist(sys); legend('Real Model Result Plot','Simulate Model Nyquist Curve'); break; end end end 例4.6 书中例 6-3 已知某对象在一组ω下的幅值和相角,试采用 Levy 拟合确定其传递函数. 频率 0.2 0.24 0.28 0.32 0.36 0.4 0.44 0.48 0.52 0.56 幅值 0.167 0.176 0.186 0.199 0.214 0.231 0.252 0.276 0.307 0.345 相角 41.4 41.16 39.81 37.67 34.98 31.87 28.38 24.47 20.06 14.97 频率 0.6 0.64 0.68 0.72 0.76 0.8 0.85 0.90 0.95 1.00 幅值 0.394 0.457 0.536 0.630 0.712 0.730 0.634 0.499 0.388 0.309 相角 8.868 1.241 -8.77 -22.36 -40.67 -62.51 -88.18 -106.7 -119 -127.4 频率 1.05 1.10 1.20 1.30 1.40 1.60 1.80 2.00 2.20 2.40 幅值 0.252 0.210 0.154 0.119 0.095 0.065 0.047 0.036 0.029 0.023 相角 -133.4 -138 -144.6 -149.2 -152.7 -157.6 -160.7 -162.7 -164 -164.8 解:首先将上表中的数据转换成实频和虚频的形式,然后调用 Levy 拟合的函数即可. w=[0.2:0.04:0.8 0.85:0.05:1.10 1.20:0.10:1.40 1.60:0.20:2.40]; A=[0.167 0.176 0.186 0.199 0.214 0.231 0.252 0.276 0.307 0.345 0.394 0.457 0.536 0.630 0.712 0.730 0.634 0.499 0.388 0.309 0.252 0.210 0.154 0.119 0.095 0.065 0.047 0.036 0.029 0.023]; f=[41.4 41.16 39.81 37.67 34.98 31.87 28.38 24.47 20.06 14.97 8.868 1.241 -8.77 -22.36 -40.67 -62.51 -88.18 -106.7 -119 -127.4 -133.4 -138 -144.6 -149.2 -152.7 -157.6 -160.7 -162.7 -164 -164.8]; %角度转化为弧度 f=f*pi/180; %转化为实频和虚频 P=A.*cos(f); Q=A.*sin(f); sigma=0.0001; %最高阶次为4 order=4; %调用levy函数 levy(w,P,Q,order,sigma); 51 得到的结果如下: -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Nyquist Diagram Real Axis Imaginary Axis Real Model Result Plot Simulate Model Nyquist Curve 二阶-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Nyquist Diagram Real Axis Imaginary Axis Real Model Result Plot Simulate Model Nyquist Curve 三阶-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Nyquist Diagram Real Axis Imaginary Axis Real Model Result Plot Simulate Model Nyquist Curve 四阶52 由上面结果可知,二阶模型拟合的结果较差,三阶、四阶模型拟合的结果非常好. 4.3 相关分析辨识系统脉冲响应函数 4.3.1 相关函数的求取当sNT T = ,t 和τ均为 s T 的整数倍时,即sskT , iT t = τ = ,i=1, 2, …, n,k=1, 2, …, n, 则自相关函数为 [ ] ∑ ∫ = ∞ → + ≈ τ + = = τ N i s s s T T s uu uu kT iT u iT u N dt t u t u T kT R R 1 0 ) ( ) ( 1 ) ( ) ( 1 lim ) ( ) ( [ ] ∑ ∫ = ∞ → + ≈ = τ + τ = τ N i s s s s uy T T uy kT iT y iT u N kT R dt t y u T R 1 0 ) ( ) ( 1 ) ( ) ( ) ( 1 lim ) ( 根据上面的公式很容易由实验数据计算出输入的自相关函数及输入输出之间的互相关函数. 需要注意的是,如果要计算的是 k=1, 2, …, n 的相关函数,则实验数据长度至少应为 n+N. 例4.7 已知某周期信号周期为 127,四个周期内的数据如下图所示,试计算其两个周期内的自相关函数. 0 100 200 300 400 500 600 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 解:根据要求 N=2*127,求i时刻自相关函数的程序段如下: Ruu(i)=0; for j=1:N Ruu(i)=Ruu(i)+u(j)*u(j+i); end Ruu(i)=Ruu(i)/N; 求得两个周期内的自相关函数如下图所示. 53 0 50 100 150 200 250 300 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 4.3.2 用伪随机信号辨识脉冲响应函数如果输入为 PRBS 信号,用相关分析法辨识系统模型的算法如下: (i) 输入试验信号的参数 a, s T , L;输入试验信号 u(t)和输出信号 y(t)的采样数据 u(k)、y(k),k=0,1,……,L-1; (ii) r=0; (iii) k=0,M=0; (iv) 计算 l=k+r; (v) 若l>L-1,转(vi)否则转(vii); (vi) l=k+r-L; (vii) M=M+ u(k)y(l); (viii) 如果 k=L-1,转(x),否则转(ix); (ix) k=k+1,转(iv); (x) L M rT R s uy = ) ( ; (xi) 如果 r=L-1,转(xiii),否则转(xii); (xii) r=r+1, 转(iii); (xiii) 计算 ∑ ? = = 1 0 ) ( L k s uy kT R W ; (xiv) 计算 [ ] W rT R T a L L rT g s uy s s + + = ) ( ) 1 ( ) ( 2 ,r=0,1,……L-1; (xv) 结束. 仿真中 PRBS 信号的生成用 idinput(),用法如下: U = IDINPUT(N,TYPE,BAND,LEVELS) 54 参数说明如下: U:返回的产生的输入信号,为列向量或 N*nu 维矩阵. N:输入信号的长度. N = [N Nu]产生 N*Nu 的输入(Nu 个输入通道); N = [P Nu M]产生 M*(P*Nu)的周期输入信号,周期为 P,M 个周期; 缺省值为 Nu = 1 和M=1. TYPE:产生的输入信号的类型,可以取下列值之一: 'RGS':随机高斯信号; 'RBS':随机二位式信号; 'PRBS':伪随机二位式信号; 'SINE':正弦和信号; Default: TYPE = 'RBS'. BAND:1*2 的行向量,定义输入频谱的频带,对'RS','RBS'和'SINE'类型 BAND = [LFR,HFR], LFR 和HFR 分别是与 Nyquist 频率的比值表示的通带的下限和上限 (通常是 0~1 之间的数) ,对'PRBS'类型,BAND = [0,B],信号在 1/B(时钟周期)间隔之间是常数.缺省值为 BAND =[0 1]. LEVELS = [MI, MA]:2*1 的行向量,定义输入的高度,对'RBS','PRBS'和'SINE', 高度调整为输入信号始终在 MI 到MA 之间,对'RGS',MI 为信号的均值减标准差,MA 为信号的均值加标准差.缺省值为 LEVELS = [-1 1]. 对'PRBS',如果 M>1,数据序列长度和周期将调整为最长 PRBS 周期的整数,如果 M=1,周期选为比 P=N 长.在多输入情况下,信号最大程度地移位.这意味着 P/Nu 是辨识此信号激励的系统的模型阶次的上界. 对'SINE',正弦从 freq=2*pi*[1:Grid_Skip:fix(P/2)]/P 和pi*[BAND(1) BAND(2)]所成的频率网格中选择(对Grid_Skip 见下面说明), 对多输入信号, 不同的输入用这个网格中的不同频率.可以通过[U,FREQS] = IDINPUT(....)返回选择的频率,FREQS 的ku 行包含第 ku 个输入的频率.产生的信号受第五个输入参数 SINEDATA 的影响: U = IDINPUT(N,TYPE,BAND,LEVELS,SINEDATA) 其中 SINEDATA= [No_of_Sinusoids, No_of_Trials, Grid_Skip],表示 No_of_Sinusoids 平均地覆盖在 BAND 上,尝试 No_of_Trials 个不同的随机的相对相位,直到找到最低幅度的信号. 缺省值为 SINEDATA=[10,10,1]. 生成输入数据后,可以利用实验三的方法进行仿真得到输出数据,在本实验中,也可以使用 matlab 中的 sim 函数.示例如下: sys = tf(1, [1 1 5]); s=tf('s'); mod = idpoly(sys); a=1; Ts=0.1; L=255; 55 u=idinput([L 1 8],'prbs',[],[-a a]); u1 = iddata([],in,Ts); y1=sim(mod,u1); 得到的 y1 为iddata 类型的数据,可以通过 y1.y 得到 double 类型的数据.根据输入输出数据利用相关分析法可以得到系统的脉冲响应函数. 例4.7 试根据上述程序生成的输入输出数据辨识对象的脉冲响应函数. 解:在上述程序中,我们产生了 8 个周期的数据,辨识对象的脉冲响应函数的过程中, 我们计算 4 个周期的互相关函数,得到的辨识结果如下图所示: 0 50 100 150 200 250 300 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 辨识结果理论计算值辨识出对象的脉冲响应函数后,可以利用 4.2.1 节的方法计算其频率特性,进而可以用 levy 法拟合出传递函数模型. 由于相关分析法有较强的抗干扰能力,如果输出信号包含噪声,也能得到较好的辨识结果,读者可在仿真中用 awgn()对输出数据添加白噪声,并通过仿真验证此结果. 实验内容 1 分别编制阶跃响应与方波响应相互转化的程序、由阶跃响应数据获得对象传递函数模型的通用程序、由阶跃响应数据获得对象微分方程模型的通用程序,并通过实例验证程序. 2 采用相关分析方法辨识下面二阶系统的脉冲响应函数, 并用 Levy 方法拟合出传递函数 25 s 5 G(s) 2 + + = s 分别考虑输出无噪声和包含输出噪声的情况 (如噪信比分别是5%和10%等情况) 进行辨识, 并对结果进行分析. 56 实验五控制系统计算机辅助分析与设计的频域方法实验目的: 掌握传递函数的频率特性计算方法及根据频域性能指标设计串联校正装置的方法. 5.1 传递函数的频率特性计算 (1) 因式分解形式的传递函数频率特性计算的算法流程如下: (i) 输入分子、分母的一阶项个数 M1、N1,二阶项个数 M2,N2;输入各个多项式系数Ti 和zi(i=1,2,…,M1) ,ai,bi 和ci(i=1,2,…,M2) ;τi 和pi(i=1,2,…,N1) ,di, ei 和fi(i=1,2,…,N2) ;输入原极点的个数 v; (ii) 输入欲计算的频率点数 n,下限频率 ωmin 和上限频率 ωmax. (iii) 令k=1, k ω = min ω , 1 min max ? ? = ? n ω ω ω ; (iv) 调用计算一阶多项式的子程序, 分别计算分子、分母一阶项的幅值 A1i(ωk)、 B1j(ωk) 和相位 Φ1i(ωk)、Ψ1j(ωk),其中 i=1,2,…,M1,j=1,2,…,N1; (v) 调用计算二阶多项式的子程序, 分别计算分子、分母二阶项的幅值 A2i(ωk)、 B2j(ωk) 和相位 Φ2i(ωk)、Ψ2j(ωk),其中 i=1,2,…,M2,j=1,2,…,N2; (vi) 计算原极点的幅值和相位: ? 90 ) ( ) ( 3 3 * = Φ = v A k v k k ω ω ω (vii) 计算总的幅值和相位;  ∏ ∏ ∏ ∏ = = = = = 2 1 2 1 1 2 1 1 1 2 1 1 ) ( ) ( ) ( ) ( ) ( N j k j N j k j v k M i k i M i k i k ω B ω B ω ω A ω A ω A ) ( ) ( ) ( ) ( ) ( ) ( 3 1 2 1 1 1 2 1 1 2 1 2 1 k N j k j N j k j M i k i M i k i k ω ω ω ω ω ω Φ ? Ψ ? Ψ ? Φ + Φ = Φ ∑ ∑ ∑ ∑ = = = = (viii) 如果 k=n,转(x); (ix) k=k+1; ω ω ω ? + = ?1 k k ;转(iv); (x) 输出总幅值 A(ωk)和总相位 ψ(ωk) (k=1,2,…,n) . 57 (2) 直接计算多项式表示的传递函数频率特性可参考上述步骤,在需要计算的频率点上根据公式(7-24)~(7-31)进行计算. 求得频率特性的幅值和相位以后,可以根据需要画出相应的 Nyquist 图和 Bode 图.在Matlab 下还有很多方法可以得到系统的频率特性,例如 nyquist()和bode()等. 对某一频率点上我们可以用下面函数进行计算: function [mg ph]=GetMagnitudeAndPhase(sys,Wc); % 输入参数:sys: 对象传递函数 % Wc: 频率 %返回参数:mg:幅值 % ph:相位 [num den]=tfdata(sys,'v'); snum=poly2sym(num); sden=poly2sym(den); sys1=snum/sden; sys1=subs(sys1,j*Wc); re=real(sys1); im=imag(sys1); % atan return value : -pi/2 - pi/2 if re>0&&im>0 ph=atan(im/re)*180/pi-360; % 1st quadrant elseif re<0&&im>0 ph=atan(im/re)*180/pi-180; % 2nd quadrant elseif re<0&&im<0 ph=atan(im/re)*180/pi-180; % 3rd quadrant else ph=atan(im/re)*180/pi; % 4th quadrant end mg=sqrt(re^2+im^2); 例5.1 试计算对象 1 ( ) 5 G s s = + 在0.1~100 ω = 上的频率特性. 解:计算各点上的幅值和相位,用semilogx()绘制图形,并与 bode()函数计算结果进行对比如下图所示: 58 10 -1 10 0 10 1 10 2 -50 -40 -30 -20 -10 10 -1 10 0 10 1 10 2 -100 -80 -60 -40 -20 0 计算值Bode()结果5.2 串联校正装置传递函数的计算对单变量最小相位系统和包含纯滞后环节的开环稳定系统, 校正装置辅助设计的主要步骤: 1 输入原系统的结构参数、设计控制系统所要求的性能指标(相稳定裕度 PM、误差系数Kv 及截止频率 ωc)以及控制系统仿真所需参数; 2 计算原系统的相位 Φc(ωc); 3 根据所要求的相位裕度 PM 计算校正装置的相位 Φc(ωc); 4 由Φc(ωc)的正负号判断采用微分或积分校正装置; 5 计算 α、zc1、pc1; 6 计算引入微分校正后系统的误差系数 K1,校正系数 Kc; 7 如果 K1

第一部分 控制系统仿真基础知识及 Matlab 在控制系统仿真中...

第一部分控制系统仿真基础知识及 Matlab 在控制系统仿真中...