Www2.ece.ohio-state.edu
Chapter 11
Matlab Simulink Simulation Testbed
In this Chapter the matlab simulink simulation test bed used is presented. The test bed can be used by students for understaning the concepts presented in this book. Instructors can use the testbed to design new project for students as homework problems.
Matlab Test Bed for Chapter 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Program to calculate Control Gains using RSP Voltage Controller and %%
%% discrete Sliding Mode Current Controller for Impulse %%
%% and to initializes matrices parameters for simulations %%
%% Author: Nanda Marwali %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
% Define fundamental output frequency;
ffun=60;
wfun=2*pi*ffun;
% Define control sampling time
Tsamp=1/60/51; %6536/20e6;
Cinv=540e-6; % Inverter capacitor filter
Linv=298e-6; % Inverter inductor filter
Cload=90e-6; % Grass capacitor
Ltrans=0.03*208*208/80e3/wfun; % 3% p.u transformer inductance
Rtrans=0; %0.03*208*208/80e3; % 3% p.u transformer resistance
tr=120/245; % Inverter To Output Turn ratio
Ilimit=3*80e3/245*sqrt(2); % 300% inverter current limit in Line-Line Amp peak
%% Define harmonics frequencies
w1=wfun;
w2=2*wfun;
w3=3*wfun;
w5=5*wfun;
w4=4*wfun;
w7=7*wfun;
T1=1/ffun;
% Define Delta-Wye Transformer Voltages and Currents Transfer Matrices
Tri_qd=3/2*tr*[1 sqrt(3) ;
-sqrt(3) 1 ];
Trv_qd=1/2*tr*[1 -sqrt(3);
sqrt(3) 1 ];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start design of discrete SM current controller
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Define the plant for the current controller
A=[zeros(2,2) eye(2,2)/(3*Cinv) ;
-1/Linv*eye(2) -0.01/Linv*eye(2)];
B=[zeros(2,2);
1/Linv*eye(2)];
E=[-Tri_qd/(3*Cinv);
zeros(2,2)];
F=[zeros(2,2);
-eye(2)];
C=[zeros(2,2) eye(2)];
D=zeros(2,2);
%% Discretize the plant for the current controller
sysc=ss(A,B,C,zeros(size(C,1),size(B,2)));
sysd=c2d(sysc,Tsamp,'zoh');
[Acurrd,Bcurrd,Ccurrd,Dcurrd]=ssdata(sysd);
CBinv=inv(Ccurrd*Bcurrd);
CA=Ccurrd*Acurrd;
CD=Ccurrd*F;
sysc=ss(A,E,C,zeros(size(C,1),size(B,2)));
sysd=c2d(sysc,Tsamp,'zoh');
[Acurrd1,Ecurrd,Ccurrd1,Dcurrd1]=ssdata(sysd);
CE=Ccurrd1*Ecurrd;
%%% End Discrete sliding mode controllers gains calculation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start design of Perfect RSP voltages controllers
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Define the true plant
Ao=[ zeros(2,2) eye(2,2)/(3*Cinv) zeros(2,2) -Tri_qd/(3*Cinv);
-eye(2,2)/Linv zeros(2,2) zeros(2,2) zeros(2,2);
zeros(2,2) zeros(2,2) zeros(2,2) eye(2,2)/Cload;
Trv_qd/Ltrans zeros(2,2) -eye(2,2)/Ltrans -Rtrans*eye(2,2)/Ltrans];
Bo=[zeros(2,2);
eye(2,2)/Linv;
zeros(2,2);
zeros(2,2)];
Co=[zeros(2,2) zeros(2,2) eye(2) zeros(2,2)];
%% Define the analog servo compensator
Ch0=zeros(2,2);
Ch1=[zeros(2,2) eye(2);
-w1^2*eye(2) zeros(2,2)];
Ch5=[zeros(2,2) eye(2);
-w5^2*eye(2) zeros(2,2)];
Ch7=[zeros(2,2) eye(2);
-w7^2*eye(2) zeros(2,2)];
Ch_star=[Ch1 zeros(size(Ch1)) zeros(size(Ch1)) ;
zeros(size(Ch1)) Ch5 zeros(size(Ch1)) ;
zeros(size(Ch1)) zeros(size(Ch1)) Ch7 ];
Bh_star=[zeros(2,2);
eye(2,2);
zeros(2,2);
eye(2,2);
zeros(2,2);
eye(2,2)];
% Discretize true plant
sysc=ss(Ao,Bo,Co,zeros(size(C,1),size(B,2)));
sysd=c2d(sysc,Tsamp,'zoh');
[Aod,Bod,Cd,Dd]=ssdata(sysd);
% Calculate equivalent plant+DSM current controller
C1=[eye(4) zeros(4,4)];
C2=[zeros(2,4) eye(2) zeros(2,2)];
Ad=Aod-Bod*(CBinv*CA*C1+CBinv*CE*C2);
Bd=Bod*CBinv;
% Discetize the controller
csysc=ss(Ch_star,Bh_star,eye(size(Ch_star,1)),zeros(size(Ch_star,1),size(Bh_star,2)));
[csysbc,T]=ssbal(csysc);
csysd=c2d(csysbc,Tsamp,'zoh');
[Acon_d,Bcon_d,Ccon_d,Dcon_d]=ssdata(csysd);
% Form the augmented equivalent plant and the servo compensator
Ad_big=[Ad zeros(size(Ad,1),size(Acon_d,2));
-Bcon_d*Cd Acon_d];
Bd_big=[Bd ; -Bcon_d*Dd];
% Define the weighting matrices
epsilon=1e-5;
Q2=eye(size(Acon_d,1));
state_W=0.2; % 80 kVA
Q1=state_W*eye(size(Ad,1));
Q=[ Q1 zeros(size(Ad,1),size(Acon_d,2));
zeros(size(Acon_d,1),size(Ad,2)) 5e5*Q2];
R=epsilon*eye(2);
% Now performed the optimal calculations of the gains
[Kd,S,E]=dlqr(Ad_big,Bd_big,Q,R);
Kd=-Kd;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create matrices parameters for the plant states equations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rfl=208*208/(80e3*0.8);
Rtrans=0.03*208*208/80e3;
Tri=[1 -2 1;
1 1 -2;
-2 1 1];
Trv=[0 0 -1;
-1 0 0;
0 -1 0];
Tri_qd0=3/2*tr*[1 sqrt(3) 0;
-sqrt(3) 1 0];
Trv_qd0=1/2*tr*[1 -sqrt(3);
sqrt(3) 1 ;
0 0 ];
Ks=2/3*[cos(0) cos(0-2*pi/3) cos(0+2*pi/3);
sin(0) sin(0-2*pi/3) sin(0+2*pi/3);
1 1 1 ];
Ksinv=inv(Ks) ;
% Define the B and D matrices
Bsim=[zeros(2,2);
eye(2,2)/Linv;
zeros(3,2);
zeros(3,2)];
Dsim=zeros(13,2);
% Full load
Rload33=[ Rfl 0 0;
0 Rfl 0;
0 0 Rfl];
Afl=[ zeros(2,2) eye(2,2)/(3*Cinv) zeros(2,3) -Tri_qd0/(3*Cinv);
-eye(2,2)/Linv zeros(2,2) zeros(2,3) zeros(2,3);
zeros(3,2) zeros(3,2) -Ks*inv(Rload33)*Ksinv/Cload eye(3,3)/Cload;
Trv_qd0/Ltrans zeros(3,2) -eye(3,3)/Ltrans -Rtrans*eye(3,3)/Ltrans];
Cfl=[ [eye(10)];
[zeros(3,4) Ks*inv(Rload33)*Ksinv zeros(3,3)]];
% No load
Rload33=[ 100 0 0;
0 100 0;
0 0 100];
Anl=[ zeros(2,2) eye(2,2)/(3*Cinv) zeros(2,3) -Tri_qd0/(3*Cinv);
-eye(2,2)/Linv zeros(2,2) zeros(2,3) zeros(2,3);
zeros(3,2) zeros(3,2) -Ks*inv(Rload33)*Ksinv/Cload eye(3,3)/Cload;
Trv_qd0/Ltrans zeros(3,2) -eye(3,3)/Ltrans -Rtrans*eye(3,3)/Ltrans];
Cnl=[ [eye(10)];
[zeros(3,4) Ks*inv(Rload33)*Ksinv zeros(3,3)]];
% Single phase loaded
Rload33=[ Rfl 0 0;
0 100 0;
0 0 100];
Ap1=[ zeros(2,2) eye(2,2)/(3*Cinv) zeros(2,3) -Tri_qd0/(3*Cinv);
-eye(2,2)/Linv zeros(2,2) zeros(2,3) zeros(2,3);
zeros(3,2) zeros(3,2) -Ks*inv(Rload33)*Ksinv/Cload eye(3,3)/Cload;
Trv_qd0/Ltrans zeros(3,2) -eye(3,3)/Ltrans -Rtrans*eye(3,3)/Ltrans];
Cp1=[ [eye(10)];
[zeros(3,4) Ks*inv(Rload33)*Ksinv zeros(3,3)]];
% Two phase loaded
Rload33=[ Rfl 0 0;
0 Rfl 0;
0 0 100];
Ap2=[ zeros(2,2) eye(2,2)/(3*Cinv) zeros(2,3) -Tri_qd0/(3*Cinv);
-eye(2,2)/Linv zeros(2,2) zeros(2,3) zeros(2,3);
zeros(3,2) zeros(3,2) -Ks*inv(Rload33)*Ksinv/Cload eye(3,3)/Cload;
Trv_qd0/Ltrans zeros(3,2) -eye(3,3)/Ltrans -Rtrans*eye(3,3)/Ltrans];
Cp2=[ [eye(10)];
[zeros(3,4) Ks*inv(Rload33)*Ksinv zeros(3,3)]];
% 500% load
Rload33=[ sqrt(0.2*Rfl) 0 0;
0 sqrt(0.2*Rfl) 0;
0 0 sqrt(0.2*Rfl)];
Aov=[ zeros(2,2) eye(2,2)/(3*Cinv) zeros(2,3) -Tri_qd0/(3*Cinv);
-eye(2,2)/Linv zeros(2,2) zeros(2,3) zeros(2,3);
zeros(3,2) zeros(3,2) -Ks*inv(Rload33)*Ksinv/Cload eye(3,3)/Cload;
Trv_qd0/Ltrans zeros(3,2) -eye(3,3)/Ltrans -Rtrans*eye(3,3)/Ltrans];
Cov=[ [eye(10)];
[zeros(3,4) Ks*inv(Rload33)*Ksinv zeros(3,3)]];
% Short circuit
Rload33=[ 0.001 0 0;
0 0.001 0;
0 0 0.001];
Asc=[ zeros(2,2) eye(2,2)/(3*Cinv) zeros(2,3) -Tri_qd0/(3*Cinv);
-eye(2,2)/Linv zeros(2,2) zeros(2,3) zeros(2,3);
zeros(3,2) zeros(3,2) -Ks*inv(Rload33)*Ksinv/Cload eye(3,3)/Cload;
Trv_qd0/Ltrans zeros(3,2) -eye(3,3)/Ltrans -Rtrans*eye(3,3)/Ltrans];
Csc=[ [eye(10)];
[zeros(3,4) Ks*inv(Rload33)*Ksinv zeros(3,3)]];
A.2 File SFUNFFT.M
function [sys, x0, str, ts] = sfunfft(t,x,u,flag,fftpts,npts,HowOften,...
offset,ts,averaging,nharmonics)
%SFUNPSD an S-function which performs spectral analysis using ffts.
% This M-file is designed to be used in a Simulink S-function block.
% It stores up a buffer of input and output points of the system
% then plots the power spectral density of the input signal.
%
% The input arguments are:
% npts: number of points to use in the fft (e.g. 128)
% HowOften: how often to plot the ffts (e.g. 64)
% offset: sample time offset (usually zeros)
% ts: how often to sample points (secs)
% averaging: whether to average the psd or not
%
% Two or three plots are given: the time history, the instantaneous psd
% the average psd.
%
% See also, FFT, SPECTRUM, SFUNTMPL, SFUNTF.
% Copyright (c) 1990-97 by The MathWorks, Inc.
% $Revision: 1.25 $
% Andrew Grace 5-30-91.
% Revised Wes Wang 4-28-93, 8-17-93.
% Revised Charlie Ko 9-26-96.
switch flag
%%%%%%%%%%%%%%%%%%
% Initialization %
%%%%%%%%%%%%%%%%%%
case 0
[sys,x0,ts] = mdlInitializeSizes(fftpts,npts,HowOften,offset,...
ts,averaging);
SetBlockCallbacks(gcbh);
%%%%%%%%%%
% Update %
%%%%%%%%%%
case 2
sys = mdlUpdate(t,x,u,fftpts,npts,HowOften,offset,ts,averaging,nharmonics);
case 3
nstates = npts + 2 + averaging * round(fftpts/2) + 1;
sys=x(nstates+1);
%%%%%%%%%%%%%%%%%%%%%%%
% GetTimeOfNextVarHit %
%%%%%%%%%%%%%%%%%%%%%%%
case 4
sys=mdlGetTimeOfNextVarHit(t,x,u,ts);
%%%%%%%%%
% Start %
%%%%%%%%%
case 'Start'
LocalBlockStartFcn
%%%%%%%%%%%%%%
% NameChange %
%%%%%%%%%%%%%%
case 'NameChange'
LocalBlockNameChangeFcn
%%%%%%%%%%%%%%%%%%%%%%%%
% CopyBlock, LoadBlock %
%%%%%%%%%%%%%%%%%%%%%%%%
case { 'CopyBlock', 'LoadBlock' }
LocalBlockLoadCopyFcn
%%%%%%%%%%%%%%%
% DeleteBlock %
%%%%%%%%%%%%%%%
case 'DeleteBlock'
LocalBlockDeleteFcn
%%%%%%%%%%%%%%%%
% DeleteFigure %
%%%%%%%%%%%%%%%%
case 'DeleteFigure'
LocalFigureDeleteFcn
%%%%%%%%%%%%%%%%
% Unused flags %
%%%%%%%%%%%%%%%%
case { 3, 9 }
sys = []; %do nothing
%%%%%%%%%%%%%%%%%%%%
% Unexpected flags %
%%%%%%%%%%%%%%%%%%%%
otherwise
if ischar(flag),
errmsg=sprintf('Unhandled flag: ''%s''', flag);
else
errmsg=sprintf('Unhandled flag: %d', flag);
end
error(errmsg);
end
% end sfunpsd
%
%============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%============================================================================
%
function [sys,x0,ts] = mdlInitializeSizes(fftpts,npts,HowOften,offset,...
ts,averaging)
nstates = npts + 2 + averaging * round(fftpts/2) + 1; % No. of discrete states
initial_states = [1; zeros(nstates-1,1)]; % Initial conditions
if HowOften > npts
error('The number of points in the buffer must exceed the plot frequency')
end
initial_states(nstates+1) = 0;
sizes = simsizes;
sizes.NumContStates = 0;
sizes.NumDiscStates = nstates+1;
sizes.NumOutputs = 1;
sizes.NumInputs = 2;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;
sys = simsizes(sizes);
x0 = initial_states;
str = [];
%ts = [ts offset];
ts=[-2 0];
% end mdlInitializeSizes
%
%=============================================================================
% mdlGetTimeOfNextVarHit
% Return the time of the next hit for this block. Note that the result is
% absolute time.
%=============================================================================
%
function sys=mdlGetTimeOfNextVarHit(t,x,u,ts)
if u(2)==0
sys = t + 1/60;
else
sys=t+ts;
end
% end mdlGetTimeOfNextVarHit
%
%===========================================================================
% mdlUpdates
% Compute Discrete State updates.
%===========================================================================
%
function [sys]=mdlUpdate(t,x,u,fftpts,npts,HowOften,offset,ts,averaging,nharmonics)
nstates = npts + 2 + averaging * round(fftpts/2) + 1; % No. of discrete states
if u(2)==0
sys=x;
return
end
%
% Find the figure handle associated with this block. If the handle
% is not valid (i.e. - the user closed the figure during simulation),
% mdlUpdates performs no operations.
%
FigHandle = GetSfunPSDFigure(gcbh);
if ~ishandle(FigHandle)
sys=[];
return
end
%
% Initialize the return argument.
%
sys = x;
%
% Increment the counter and store the current input in the
% discrete state referenced by this counter.
%
x(1) = x(1) + 1;
sys(x(1)) = u(1);
if rem(x(1),HowOften)==0
ud = get(FigHandle,'UserData');
buffer = [sys(x(1)+1:npts+1);sys(2:x(1))];
n = fftpts/2;
freq =(1/ts); % Multiply by 2*pi to get radians
w = freq*(0:n-1)./fftpts;
% Detrend the data: remove best straight line fit
a = [(1:npts)'/npts ones(npts,1)];
y = buffer;% - a*(a\buffer);
% Hammining window to remove transient effects at the
% beginning and end of the time sequence.
nw = min(fftpts, npts);
win = .5*(1 - cos(2*pi*(1:nw)'/(nw+1)));
g = fft(y(1:nw));
% g(1)=0.0;
% Perform averaging with overlap if number of fftpts
% is less than buffer
% ng = fftpts;
% while (ng ................
................
In order to avoid copyright disputes, this page is only a partial summary.
To fulfill the demand for quickly locating and searching documents.
It is intelligent file search solution for home and business.