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.

Google Online Preview   Download