%--------H无穷控制1.  主文件： design.m

clear                     %--清除存在内存里的数据

format short e            %----MATLAB可以将计算结果以不同的精度输出.有效数字16位加3位指数

w=logspace(-2,6,100);    %用法：logspace(a,b,n)，其中a、b、n分别表示开始值（10^a）、结束值(10^b)、元素个数(n)。

                        %功能：生成从10的a次方到10的b次方之间按对数等分的n个元素的行向量。n如果省略，则默认值为50。

%----设一般被控对象

gplant

%----设计控制器

disp(['Hoo controller Design with hinfsym'])      %matlab中disp()就是屏幕输出函数，类似于c语言中的printf（）函数

K=[]

controller;

disp('Hit any key')

pause                                     %一般是为了动态观察变化过程 pause（a）暂停a秒后执行下一条指令

if isempty(K)==1                          %isempty(msg)表示判断msg变量是否是空，或者是0。如果是真，就返回1，否则返回0。

   return

end

%----显示控制器的Bode图

contplot;

disp('Hit any key')

pause

%----仿真

simu

disp('Hit any key')

pause

%close all

%clc

%-------2.一般被控对象的设定子程序： gplant.m

Psys=nd2sys([1],[1 0 0],3.87e7);   %function sys = nd2sys(num,den,gain)

                                 %  ND2SYS converts a numerator/denominator single-input/

                                   %single-output transfer function into a SYSTEM matrix.

                                   %GAIN is an optional argument that scales the output SYS.

                                   %将某某写成系统矩阵的形式

W1sys=nd2sys([1 125.7],[1 1e-4],5e-1);

W2num1=[1 1e+4 5.7e+7];

W2dencof=1;

W2den1=[1 1.2e+4 4.04e+8];

W2sys1=nd2sys(W2num1,W2den1);

W2sys=mmult(W2sys1,W2sys1,23.9);   %function out = mmult(mat1,mat2,mat3,...,matN)

                                   %Multiplies SYSTEM/VARYING/CONSTANT matrices, currently

                                  % limited to 9 input arguments

                                  %out = mat1*mat2*mat3* *** *matN

                                   %G的输入

% W2sys=mmult(W2sys1,23.9);

W3sys=0.1;

W4sys=nd2sys([10 10*0.5*0.5e4],[1 10*5e4]);

%-----显示加权函数的Bode图

Pfr=frsp(Psys,w);     %function out = frsp(sys,omega,T,balflg)

                    % Complex frequency response of a SYSTEM matrix, produces

                     %a VARYING matrix which contains its frequency response:

                    % SYS     - SYSTEM matrix

                   % OMEGA   - response calculated at these frequencies, can

                     %also be another VARYING matrix, and then these

                 % independent variables are used

               % T       - 0 (default) continuous system, if T ~= 0 then

              % discrete evaluation with sample time T is used.

             %  BALFLG  - 0 (default) balances A matrix prior to eval, nonzero

             %    value for BALFLG leaves state-space unchanged

             % OUT     - VARYING frequency response matrix

W1fr=frsp(W1sys,w);

W2fr=frsp(W2sys,w);

W3fr=frsp(W3sys,w);

W4fr=frsp(W4sys,w);

figure(1);

vplot('liv,lm',W1fr,W2fr,W3fr,W4fr);    %function vplot([plot_type],vmat1,vmat2,vmat3, ...)Plot one or more varying matrices.  The syntax is

                                        % the same as the MATLAB plot command except that all data

                                        %is contained in VMATi, and the axes are specified by PLOT_TYPE.

                                        %The (optional) plot_type argument must be one of:'liv,lm'     log(magnitude) .vs. log(independent variable)

title('Weightign Functions');

xlabel('Frequency[rad/s]')

ylabel('Gain')

grid;

%-----设一般被控对象

systemnames='Psys W1sys W2sys W3sys W4sys';

inputvar='[w1;w2;u]';

outputvar='[W1sys;W2sys;W4sys;w2+Psys]';

input_to_W3sys='[w1]';

input_to_Psys='[u+W3sys]';

input_to_W1sys='[Psys+w2]';

input_to_W2sys='[Psys]';

input_to_W4sys='[u]';

sysoutname='Gpsys';

sysic;

%-----定义控制对象和加权函数

%------------3.控制器的设计子程序： controller.m

glow=0;ghigh=10;tol=1e-2;

while isempty(K)==1 & (ghigh<1e6);

   [K,CL,gopt]=hinfsyn(Gpsys,1,1,glow,ghigh,tol);     %见鲁棒迭代笔记55

   if isempty(K)==1 & ghigh<1e6;

       ghigh=ghigh*10;

       tol=ghigh*1e-3;

   end

   clear CL ghigh glow tol

end;%

K

%------------4.控制器的Bode图的显示子程序：contplot.m

Kfr=frsp(K,w);

figure(2);

subplot(2,1,1)      %使用方法：subplot（m,n,p）或者subplot（m n p）。

                   %subplot是将多个图画到一个平面上的工具。其中，m表示是图排成m行，n表示图排成n列，也就是整个figure中有n个图是排成一行的，一共m行，如果m=2就是表示2行图。

                   %p表示图所在的位置，p=1表示从左到右从上到下的第一个位置。

vplot('liv,lm',Kfr);

title('Bode plot of Controller')

xlabel('Frequency[rad/s]')

ylabel('Gain')

grid;

subplot(2,1,2)

vplot('liv,p',Kfr);

xlabel('Frequency[rad/s]')

ylabel('degrees')

grid;

%==============================================================

%-------------5.闭环系统的设计子程序：closedlp.m

systemnames='Psys K';

inputvar='[w;r]';

outputvar='[Psys+r;K]';

input_to_Psys='[K+w]';

input_to_K='[Psys+r]';

sysoutname='Clsys';

cleanupsysic='yes';

sysic;

%----------------6.仿真子程序：simu.m

%----仿真（H无穷控制器）

closedlp;

[pa,pb,pc,pd]=unpck(Psys);      %%直接提取出系统的状态方程模型

[ka,kb,kc,kd]=unpck(K);

[da,db,dc,dd]=unpck(sel(Clsys,1,1));

[ua,ub,uc,ud]=unpck(sel(Clsys,2,1));

%----干扰响应

t=[0:0.0001:0.1]

[y,x,t]=step(da,db,dc,dd,1,t); %求取系统单位阶跃响应：step()

                                  %y=step(num,den,t)：其中num和den分别为系统传递函数描述中的分子和分母多项式系数，t为选定的仿真时间向量，

                                  %一般可以由t=0:step:end等步长地产生出来。该函数返回值y为系统在仿真时刻各个输出所组成的矩阵

                                  %[y,x,t]=step(num,den)：此时时间向量t由系统模型的特性自动生成, 状态变量x返回为空矩阵。

                                   %[y,x,t]=step(A,B,C,D,iu)：其中A,B,C,D为系统的状态空间描述矩阵，iu用来指明输入变量的序号。x为系统返回的状态轨迹。

[u,x,t]=step(ua,ub,uc,ud,1,t);

figure(3);

subplot(2,1,1)

plot(t,y)

axis([0 0.1 0 16])

xlabel('Time[s]')

ylabel('Amplitude')

title('Step disturbance response')

grid on;

%---控制输入

subplot(2,1,2)

figure(4)

plot(t,u)

axis([0 0.01 -1.7 0])

xlabel('Time[s]')

ylabel('u[V]')

title('Input')

grid on