pmc.m

MATLAB预测控制程序实例

%%%%%%%%%%%%%%%%%%%
%%Describe System%%
%%%%%%%%%%%%%%%%%%%
%DC-motor with elastic shaft
%Parameters(MKS)
Lshaft=1.0;                
dshaft=0.02;                
shaftho=7850;               
G=81500*1e6;               
tauam=50*1e6;               
Mmotor=100;               
Rmotor=0.1;               
Jmotor=0.5*Mmotor*Rmotor^2;       
Bmotor=0.1;
R=20;
Kt=10;
gear=20;
Jload=50*Jmotor;
Bload=25;
Ip=pi/32*dshaft^4;
%(circular)section
Kth=G*Ip/Lshaft;
%(Torque/angle)
Vshaft=pi*(dshaft^2)/4*Lshaft;
Mshaft=shaftho*Vshaft;
Jshaft=Mshaft*0.5*(dshaft^2/4);
JM=Jmotor;
JL=Jload+Jshaft;
Vmax=tauam*pi*dshaft^3/16;
Vmin=-Vmax;

%Input/State/Output
AA=[0             1                0                  0;...
    -Kth/JL    -Bload/JL       Kth/(gear*JL)          0;... 
    0             0                0                  1;...
    Kth/(JM*gear) 0        -Kth/(JM*gear^2)   -(Bmotor+Kt^2/R)/JM];
BB=[0;0;0;Kt/(R*JM)];
Hyd=[1  0  0  0];
Hvd=[Kth 0 -Kth/gear 0];
Dyd=0;
Dvd=0;
sys=ss(AA,BB,[Hyd;Hvd],[Dyd;Dvd]);  %得到被控系统LTI离散时间模型

%%%%%%%%%%%%%%%%%%%%%
%%Controller Desgin%%
%%%%%%%%%%%%%%%%%%%%%
ManipulatedVariables=struct('Min',-220,'Max',220,'Units','V');  %设置输入限制
OutputVariables(1)=struct('Min',-Inf,'Max',Inf,'Units','rad');  %设置输出限制1
OutputVariables(2)=struct('Min',-78.5,'Max',78.5,'Units','Nm'); %设置输出限制2
Weights=struct('Input',0,'InputRate',0.05,'Output',[10 0]);     %设置权系数
Modle.Plant=sys;
Modle.Plant.OutputGroup={[1],'Measured';[2],'Unmeasured'};
Ts=0.1;                                                    %设置采样时间
PredictionHorizon=12;                                      %设置预测时域
ControlHorizon=2;                                          %设置控制时域

%%%%%%%%%%%%%%
%%MPC Plant%%%
%%%%%%%%%%%%%%
ServoMPC=mpc(Modle,Ts,PredictionHorizon,ControlHorizon,Weights,ManipulatedVariables,OutputVariables);
%get(ServoMPC)                                     %查看控制器属性
%display(ServoMPC) %get plant attribute            %查看控制器属性
%get(ServoMPC,'OutputVariables')                   %查看控制器参数
%set(ServoMPC,'PredictionHorizon'，10)             %修改控制器参数

%%%%%%%%%%%%%%
%%Simulation%%
%%%%%%%%%%%%%%
TimeSteps=round(10/Ts);
r=[pi  0];                           %输入角位置参考Angle=pi
[y,t,u]=sim(ServoMPC,TimeSteps,r);
figure(1)
%subplot(311)
plot(t,y(:,1),[0  t(end)],pi*[1  1])
title('Angular Position(radians)')
xlabel('Elapsed Time(seconds)')
grid on
figure(2)
%subplot(312)
plot(t,y(:,2),[0  t(end)],[-78.5  -78.5])
title('Torque(nM)')
xlabel('Elapsed Time(seconds)')
grid on
figure(3)
%subplot(313)
stairs(t,u)
title('Applied Voltage(volts)')
xlabel('Elapsed Time(seconds)')
grid on
figure(4)
subplot(311)
plot(t,y(:,1),[0  t(end)],pi*[1  1])
title('Angular Position(radians)')
xlabel('Elapsed Time(seconds)')
grid on
subplot(312)
plot(t,y(:,2),[0  t(end)],[-78.5  -78.5])
title('Torque(nM)')
xlabel('Elapsed Time(seconds)')
grid on
subplot(313)
stairs(t,u)
title('Applied Voltage(volts)')
xlabel('Elapsed Time(seconds)')
grid on
% U=[t';u']'
% Y=[t';y']'