- 积分
- 32
- 注册时间
- 2009-7-29
- 仿真币
-
- 最后登录
- 1970-1-1
|
发表于 2010-11-24 15:15:30
|
显示全部楼层
来自 湖南湘潭
1、可以采用函数参数传递的方式,避免每执行一次子函数就重新赋值一次(重复无用),如P,C,H,Q,R。
TO:lin2009,初始化那个函数里是可以自己加参数的么
function [sys,x0,str,ts,P]=mdlInitializeSizes
s-functon函数及子函数的输入、输出都有固定的格式。
输出不能改变,如初始化子函数的输出就是 [sys,x0,str,ts],不能改变,不能自己加参数。
但是输入可以自己加参数。
因此,你的一些参数如P,C,H,Q,R等与仿真系统变量无关的常数项,可以考虑直接利用参数传递进来。如function [sys,x0,str,ts] = EKFLITER(t,x,u,flag)
可以改为function [sys,x0,str,ts] = EKFLITER(t,x,u,flag,P,C,H,Q,R)
2、若不通过上面用函数参数传递的方式。在各子函数之间共用一个变量(即公共变量),则需要声明是gobal类型或persistent类型。
参考程序:- function [ sys, x0, str, ts ] = EKFLITER(t, x, u, flag)
- %This S - Function is to indentify the inertia online using extended kalman
- %filter method.
- %x(1) = w(estimate speed)
- %x(2) = TL(load torque)
- %x(3) = 1/J(J is inertia of motor)
- %u(1) = Te
- %u(2) = w(the real speed)
- switch flag,
- case 0,
- [ sys, x0, str, ts ] = mdlInitializeSizes;
- case 2,
- sys = mdlUpdate(t, x, u);
- case 3,
- sys = mdlOutputs(t, x, u);
- case {1, 4, 9},
- sys = [ ];
- otherwise
- error([ 'Unhandled flag = ', num2str(flag) ]);
- end
- function [ sys, x0, str, ts ] = mdlInitializeSizes
- sizes = simsizes;
- sizes.NumContStates = 0;
- sizes.NumDiscStates = 3;
- sizes.NumOutputs = 2;
- sizes.NumInputs = 2;
- sizes.DirFeedthrough = 0;
- sizes.NumSampleTimes = 1;
- sys = simsizes(sizes);
- x0 = [ 0 0 1500 ];
- str = [ ];
- ts = 1e-4;
- function sys = mdlUpdate(t, x, u)
- ts = 1e-4;
- A = [ 1 - x(3)*ts 0; 0 1 0; 0 0 1 ];
- B = [ x(3)*ts 0 0 ]';
- % P = diag([ 1 1 10 ]);
- C = [ 1 0 0 ];
- H = [ 1 0 0 ];
- Q = diag([ 1e-5 1e-4 5e-6 ]);
- R = 5e-3;
- y = u(2);
- Y = C*x;
- G = [ 1 - x(3)*ts u(1) - x(2)*ts; 0 1 0; 0 0 1 ];
- %The first step
- x = A*x + B*u(1);
- persistent P
- if isempty(P)
- P = diag([ 1 1 10 ]);
- else
- P = G*P*G' + Q;
- end
- % 可以将C H Q R 等变量按上面的方式定义为persistent。
- % 注意顺序
- %
- % P = G*P*G' + Q;%这里应该是用P = G*Pe*G’ + Q但是由于Pe在之前没有定义所以不能用, 现在的问题在于, 本来是一个迭代循环, 可是现在P每次更新都被重置为diag[ 1 1 10 ], 所以程序的收敛速度非常慢。看看哪位朋友能帮忙解决一下, 就是起始的P值为diag[ 1 1 10 ], 之后P = G*Pe*G’ + Q, 多谢大家了, 已经弄了好几天了, 没有找到方法。
- % %The second step
- K = P*H'/(H*P*H' + R);
- Xe = x + K*( y - Y );
- Pe = P - K*H*P;
- %Result
- sys(1) = Xe(1);
- sys(2) = Xe(2);
- sys(3) = Xe(3);
- function sys = mdlOutputs(t, Xe, u)
- %Output sys(1) is load torque
- %sys(2) is inertia J
- sys(1) = Xe(2);
- sys(2) = 1/Xe(3);
复制代码 |
评分
-
1
查看全部评分
-
|