H_infinity_Design.rar_H 控制器_H-infinity_H无穷鲁棒控制_matlab H控制器_鲁棒H∞控

%__________________________________________________________ % % tank_DKiteration.m % % Perform an H_infinity/mu synthesis design for a process % tank system. The main purpose of this code is to % demonstrate the use of the software for robust H_infinity % control design. % % Requires version 3.0 (Sept/2004) of the robust control % toolbox. % % Roy Smith, 23/May/2006 % %__________________________________________________________ % The system consists of a water tank, with controlled hot and % cold water inlet flows. An outlet flow is located at the % bottom of the tank. The measured variables of interest % are the level and temperature of the water in the tank. % The tank system, and the perturbation model is described in: % % "Model validation for robust control: an experimental process % control application," Roy S. Smith, Automatica, Vol. 31, % No. 11, pp. 1637-1647, 1995. % _________________________________________________________ clear all % Nominal model: % Physical constants: (scaled dimensionless units) A1 = 91.4; % tankcross-sectional area alpha = 1.34; % exit flow gain beta = 0.6; % exit flow bias th = 1.0; % hot water temperature tc = 0.0; % cold water temperature act_tc = 0.05; % actuator time constant % Nominal operating point: h1nom = 0.47; t1nom = 0.5; % Linearized model with % input variables: [fh; fc] % state variables: [E; f1] % output variables: [h1; t1] A_tank = [-(1+beta/h1nom)/(A1*alpha), beta*t1nom/(A1*h1nom); 0, -1/(A1*alpha) ]; B_tank = [th/A1, tc/A1; 1/(A1*alpha), 1/(A1*alpha)]; C_tank = [0, alpha; 1/h1nom, -t1nom*alpha/h1nom]; D_tank = [0, 0; 0, 0]; P_tank = ss(A_tank,B_tank,C_tank,D_tank); % Actuator model. % We use a model of the form: exp(-theta*s)*k/(tau*s + 1) % with 0.8 <= k <= 1.2, 1.5 <= tau <= 2.5, and 0.5 <= theta <= 1.0. % A second order Pade approximation models the delay. s = tf('s'); P_act = (s^2 - 8*s + 21.3)/((2*s+1)*(s^2+8*s+21.3)); % Examine the frequency response. omega = logspace(-4,2,250); P_tank_f = frd(P_tank,omega); P_act_f = frd(P_act,omega); figure(1) subplot(2,1,1) loglog(abs(P_tank_f(1,1)),abs(P_tank_f(1,2)),... abs(P_tank_f(2,1)),abs(P_tank_f(2,2)),... abs(P_act_f)); axis([min(omega),max(omega),1e-5,10]) grid legend('h1 <- fh','h1 <- fc','t1 <- fh','t1 <- fc','P_{act}') xlabel('Frequency [rad/sec]') ylabel('Magnitude') r2d = 180/pi; subplot(2,1,2) semilogx(r2d*unwrap(angle(P_tank_f(1,1))),r2d*unwrap(angle(P_tank_f(1,2))),... r2d*unwrap(angle(P_tank_f(2,1))),r2d*unwrap(angle(P_tank_f(2,2))),... r2d*unwrap(angle(P_act_f))) grid xlabel('Frequency [rad/sec]') ylabel('Phase [deg]') title('Nominal model') % -------------- Perturbation weights --------------------- % Actuator model (see Laughlin et al., 1987). W_act = 0.35*(3.83*s^2 + 9.4*s + 1.6)/(1.5*s^2 + 13*s + 8); % Output multiplicative perturbation weights W_h1m = 0.5*s/(1 + 0.25*s); % height response W_t1m = (0.1 + 16.5*s)/(1+0.2*s); % temperature response % -------------- performance weights ---------------------- % Noise weights W_hn = 0.01; W_tn = 0.03; W_hn = ss(W_hn); % convert to state-space W_tn = ss(W_tn); % Disturbance weights W_h1d = 2.5/(500*s+1); % low frequency disturbance W_t1d = 5/(500*s+1); % low frequency disturbance % Performance weights W_h1perf = 20/(600*s+1); W_t1perf = 20/(800*s+1); % Actuator penalties W_fhperf = 0.04*(1+40*s)/(1+0.1*s); W_fcperf = 0.04*(1+40*s)/(1+0.1*s); % Frequency response of the weights W_h1m_f = frd(W_h1m,omega); W_t1m_f = frd(W_t1m,omega); W_act_f = frd(W_act,omega); W_hn_f = frd(W_hn,omega); W_tn_f = frd(W_tn,omega); W_h1d_f = frd(W_h1d,omega); W_t1d_f = frd(W_t1d,omega); W_h1perf_f = frd(W_h1perf,omega); W_t1perf_f = frd(W_t1perf,omega); W_fhperf_f = frd(W_fhperf,omega); W_fcperf_f = frd(W_fcperf,omega); figure(2) subplot(1,1,1) loglog(abs(W_h1m_f),'r-',abs(W_t1m_f),'m-',... abs(W_act_f),'r--',... abs(W_hn_f),'g--',abs(W_tn_f),'b--',... abs(W_h1d_f),'c-',abs(W_t1d_f),'y-',... abs(W_h1perf_f),'g-',abs(W_t1perf_f),'b-',... abs(W_fhperf_f),'c-',abs(W_fcperf_f),'c-.') grid legend('W_{h1m}','W_{t1m}','W_{act}','W_{hn}','W_{tn}',... 'W_{h1d}','W_{t1d}',... 'W_{h1perf}','W_{t1perf}','W_{fhperf}','W_{fcperf}') xlabel('Frequency [rad/sec]') ylabel('Magnitude') title('Perturbation weights') % ---------------- interconnection structure ------------ % The interconnection structure has the following i/o % connections. The numbers corrrespond to the port % numbers on the simulink block diagram. The letters % correspond to: % ________ % | | % Delta inputs: z <---| |<--- v: Delta output % | | % cost/errors: e <---| P |<--- w: exogenous inputs % | | % measurements: y <---| |<--- u: control actuation % |_______| % % % Outputs Inputs % % Delta1 (hot act) z1 v1 % Delta2 (cold act) z2 v2 % Delta3 (h1 pert) z3 v3 % Delta4 (t1 pert) z4 v4 % % h1 tracking err e5 w5 h1 noise % t1 tracking err e6 w6 t1 noise % fh act penalty e7 w7 h1 reference % fc act penalty e8 w8 t1 reference % w9 h1 disturbance % w10 t1 disturbance % % h1 reference y9 u11 fh command % t1 reference y10 u12 fc command % h1 measurement y11 % t1 measurement y12 % % The controller will have 4 inputs and 2 outputs. Note % the positive feedback is assumed in the implementation % of the controller. [A_P,B_P,C_P,D_P] = linmod('tank_model'); P = ss(A_P,B_P,C_P,D_P); % Create indices for each block. Iz = [1:4]'; Iv = [1:4]'; Ie = [5:8]'; Iw = [5:10]'; Iy = [9:12]'; Iu = [11:12]'; % define dimensions for each nz = length(Iz); nv = length(Iv); ne = length(Ie); nw = length(Iw); nmeas = length(Iy); nctrl = length(Iu); % Define block structures for each problem. See the % mussv documentation for the format. RS_blk = [1,1; 1,1; 1,1; 1,1]; NP_blk = [6,4]; RP_blk = [RS_blk;NP_blk]; % ---------------- nominal design ----------------------- Pnomdesign = P([Ie;Iy],[Iw;Iu]); % select [e;y] <- [w;u] Pnomdesign = minreal(Pnomdesign); % remove unobservable/uncontrollable % states. [Knom,Gnom,gamma,info] = hinfsyn(Pnomdesign,nmeas,nctrl,... 'METHOD','ric',... % Riccati solution 'DISPLAY','on',... % verbose 'TOLGAM',0.1); % gamma tolerance % Check this design and see if we are happy with this as the best % possible nominal performance. At this stage we should retune % the weights depending of frequency responses, step response, % etc. clp_eig = eig(Gnom); fprintf('Max real part weighted closed loop eig = %g',... max(real(clp_eig))); if max(real(clp_eig)) > -eps, fprintf(' UNSTABLE !\n\n'); else fprintf(' (stable)\n\n'); end % Set up an unweighted simulation to examine this. K = Knom; % set controller for simulink block [Aclp,Bclp,Cclp,Dclp] = linmod('tank_unweighted_clp'); Gclp = ss(Aclp,Bclp,Cclp,Dclp); % compare the nominal closed-loop transfer functions Gclp_nom = Gclp([5:8],[5:10]); Gclp_nom = minreal(Gclp_no