1
Modelowanie w Matlabie
Program symulacyjny  jego funkcje i struktura
Komputerową symulację polega na zrealizowaniu określonych obliczeń przez specjalnie do tego
celu przeznaczony program. Mimo, że konkretna symulacja, w tym symulacja układów
energoelektronicznych, może dotyczyć bardzo różnorodnych przypadków, to jednak można mówić o
grupie pewnych, typowych funkcji, jakie taki program zazwyczaj realizuje.
Funkcje programu symulacyjnego:
" inicjacja i realizacja procesu obliczeniowego przebiegów czasowych modelu,
" inicjacja i realizacja procesu archiwizacji rezultatów symulacji,
" prezentacja rezultatów symulacji,
" funkcje opcjonalne.
Inicjacja procesu obliczeniowego przebiegów symulowanego układu, polega na
zadeklarowaniu wymaganych przez niego danych i wykonanie odpowiednich instrukcji
inicjacyjnych. W przypadkach korzystających z typowego solwera układu zwyczajnych
równań różniczkowych sprowadza się to zainicjowanie jego działania. W Matlabie dostępne
solwery sÄ… funkcyjnymi m-plikami o nazwach rozpoczynajÄ…cych siÄ™ literami ode (ode45,
ode23, ode113, ode15s, ode23s, ode23t, or ode23tb ) a będącymi skrótem w j. angielskim
(ordinary differential equations) oznaczającym zwyczajne równania różniczkowe.
Podstawowa forma odwołania do wszystkie tych solwerów jest taka sama przyjmując
przyjmując postać jednej z następujących instrukcji Matlaba
Syntaktyka solvera (podstawowa):
[T,Y] = solver(odefun,tspan,y0)
[T,Y] = solver(odefun,tspan,y0,options)
[T,Y] = solver(odefun,tspan,y0,options,p1,p2,,..)
gdzie:
solver: ode45, ode23, ode113, ode15s, ode23s, ode23t, or ode23tb
odefun  uchwyt funkcji, która wywoływana zwraca wektor pochodnych dy/dt
np. przy wywołaniu:
dy/dt = mod(t,y)
lub
dy/dt = mod(t,y, p1,p2,& )
oraz zastÄ…pieniu nazwy odefun przez nazwÄ™ 'model' przez instrukcjÄ™
model = @mod
określa się uchwyt funkcji model.
 2
np.
[T,Y] = od23('mod',tspan,y0,options,p1,p2,,..)
[T,Y] = od23(model,tspan,y0,options,p1,p2,,..)
tspan  wektor chwil czasu,
dla tspan = [tspan(1; tspan(2)]  przedział czasu rozwiązania (wektora T)
dla tspan = [tspan(1); tspan(2); & ; tspan(n)]  wektor chwil czasu równy
wektorowi czasu rozwiÄ…zania (wektora T)
p1,p2,,..  parametry
tpan(1), y0  warunek poczÄ…tkowy
[T,Y]  numeryczne rozwiÄ…zanie
Syntaktyka ustawiania opcji solvera:
options = odeset('name1',value1,'name2',value2,...)
Wybrane ważniejsze opcje (uzyskane przez zrealizowanie instrukcji odeset):
odeset
AbsTol: [ positive scalar or vector {1e-6} ]
RelTol: [ positive scalar {1e-3} ]
OutputFcn: [ function_handle ]
Refine: [ positive integer ]
InitialStep: [ positive scalar ]
MaxStep: [ positive scalar ]
odeset
AbsTol: [ positive scalar or vector {1e-6} ]
RelTol: [ positive scalar {1e-3} ]
NormControl: [ on | {off} ]
NonNegative: [ vector of integers ]
OutputFcn: [ function_handle ]
OutputSel: [ vector of integers ]
Refine: [ positive integer ]
Stats: [ on | {off} ]
 3
InitialStep: [ positive scalar ]
MaxStep: [ positive scalar ]
BDF: [ on | {off} ]
MaxOrder: [ 1 | 2 | 3 | 4 | {5} ]
Jacobian: [ matrix | function_handle ]
JPattern: [ sparse matrix ]
Vectorized: [ on | {off} ]
Mass: [ matrix | function_handle ]
MStateDependence: [ none | {weak} | strong ]
MvPattern: [ sparse matrix ]
MassSingular: [ yes | no | {maybe} ]
InitialSlope: [ vector ]
Events: [ function_handle ]
Syntaktyka odczytu aktualnych, (wskazanych) opcji solvera:
odeget(options,'name1',,'name2',...)
Przykład:
options = odeset('AbsTol',1e-7, 'RelTol',1e-3, 'OutputFcn',@wyfun, 'MaxStep',hmx);
%options = odeset('AbsTol',1e-6, 'RelTol',1e-3, 'OutputFcn',@wyfun);% 'MaxStep',hmx);
% r - rzÄ…d modelu (rozmiar wektora stanu)
model = @mod; r = 2; % Function Handle dyskretnego modelu z komutacjÄ…
x0 = zeros(1,r); % - wektor wierszowy
p2 = model; % dane o modelu dla 'OutputFun'
p3 = []; % parametr zarezerwowany dla 'OutputFcn'
p4 = tp; % czas poczÄ…tkowy tp (dla ustawiania stanu poczÄ…tkowego zmiennych
logicznych)
p5 =[]; % rezerwa
[t,x] = ode23(model,tspan,x0,options,p1,p2,p3,p4,p5);
 4
Przykład programu PS1  Model układu regulacji prądu wzbudzenia
Iwz eIw Regulator ust Przeksz- u Iw
P
Obwód
P/PI tałtnik RwLw
Tor
pomiaru
prÄ…du
Rys. 1 Schemat blokowy układu regulacji prądu wzbudzenia
Człon proporcjonalno-całkujący z ograniczeniami
Rp
iwe Rwe
i-
uwy
uwe
Cs Rs is
Rp uCs uRs
iogr
u-~= 0
umx umi
uogr
Analogowy człon PI z ograniczeniami
 5
uwy( t )= "{ }
= "{ }
= , "{ }
Dla = uogr uogr"{umin ,umx} (22)
iogr( t ) `" 0 (23)
`"
`"
`"
StÄ…d
uogr = uRs( t ) + uCs( t ) (24)
= +
= +
= +
duCs( t )
uogr = Rsis( t ) + uCs( t ) = RsCs + uCs( t ) (25)
= + = +
= + = +
= + = +
dt
duCs ( t ) 1
= (uogr - uCs ( t ))
= ( - )
= ( - )
= ( - )
(26)
dt RsCs
Rp
Rp
Rwe
iwe Rwe i- uwy
-uwy
Rwe
uwe
Cs Rs is
Rp
Rp uCs uRs
iogr
u-~= 0
umx umin
uogr
Rys.2b Analogowy człon PI z ograniczeniami i inwertorem sygnału wyjściowego
Jeśli wprowadzi się inwertor sygnału wyjściowego, co odpowiada poniższym
podstawieniom:
= - = - = - = -
= - = - = - = -
= - = - = - = -
Ć Ć Ć Ć
uwy( t ) = -uwy( t ), uCs( t ) = -uc( t ), umx = -umin , umin = -umx (27)
to otrzyma się następujący opis matematyczny członu PI z ograniczeniami
 6
umx K uwe( t ) + uc ( t ) e" umx
+ e"
Å„Å‚ + e"
Å„Å‚ + e"
Å„Å‚
Å„Å‚
R
ôÅ‚K
ôÅ‚
ôÅ‚
ôÅ‚
uwy ( t ) = uwe ( t ) + uc ( t ) K uwe ( t ) + uc( t )" ( umin ,umx ) (28a)
= + + "
= + + "
= + + "
òÅ‚
òÅ‚
òÅ‚
òÅ‚
R R
ôÅ‚u
ôÅ‚
ôÅ‚
ôÅ‚
KRuwe ( t ) + uc( t ) d" umin
+ d"
+ d"
+ d"
ół min
ół
ół
ół
ci
Å„Å‚
Å„Å‚
Å„Å‚
Å„Å‚
(umx - uc( t )) KRuwe( t ) + uc( t ) e" umx
( - ) + e"
( - ) + e"
( - ) + e"
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
Kr
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
duc( t ) ôÅ‚c
ôÅ‚
ôÅ‚
ôÅ‚
= uwe( t ) KRuwe( t ) + uc( t )" umin ,umx ) (28b)
= + "(
= + "
= + "
òÅ‚
òÅ‚
òÅ‚
òÅ‚
i
dt
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
ci
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚
(umin - uc( t )) KRuwe( t ) + uc( t ) d" umin
( - ) + d"
( - ) + d"
( - ) + d"
ôÅ‚ Kr
ôÅ‚
ôÅ‚
ôÅ‚
ół
ół
ół
ół
Rs 1 1
= = =
gdzie KR = , ci = =
= = =
= = =
Rwe RweCws Ti
Przekształtnik:
=
gdzie - współczynnik wzmocnienia
Obwód wzbudzenia:
 bez uwzględniania nieliniowości przewodzenia przekształtnika
= -
 z uwzględnieniem nieliniowości przewodzenia przekształtnika (wersja 1)
- , e"

=
, < 0
 z uwzględnieniem nieliniowości przewodzenia przekształtnika (wersja 2)
- , e" (" - > 0

=
, < 0 (" - d"
 7
Przyklad programu PS1b  Model układu regulacji prądu twornika
z przekształtnikiem nienawrotnym
Itz Człon Itwz eIt Regulator ust Przeksz- UPt It
Obwód
tałtnik
zadajÄ…cy P/PI RL
nienawr.
-E
Tor
pomiaru
prÄ…du
Człon zadający wzorzec prędkości (prądu) z ograniczeniem 1-szej pochodnej jego
zmian
avzmin avzmx
vz
awz
vwz
Kr t
Rys. 4 Schemat blokowy członu formującego wzorzec prędkości z ograniczeniem przyspieszenia
d
(vwz ( t ))= awz ( t ) (1)
( )=
( )=
( )=
dt
 8
avz min ( - )d"
( - )d"
Å„Å‚ Kr(vwz ( t )- vz ( t ))d" avz min
Å„Å‚ ( - )d"
Å„Å‚
Å„Å‚
ôÅ‚
ôÅ‚
ôÅ‚K
ôÅ‚
avz ( t ) = (vwz ( t ) - vz ( t )) avz min d" Kr (vwz ( t )- vz ( t ))d" avzmx (2)
= ( - ) d" ( - )d"
= ( - ) d" ( - )d"
= ( - ) d" ( - )d"
òÅ‚
òÅ‚
òÅ‚
òÅ‚
r
ôÅ‚a
ôÅ‚
ôÅ‚
ôÅ‚
Kr(vwz ( t )- vz ( t ))e" avzmx
( - )e"
( - )e"
( - )e"
vzmx
ół
ół
ół
ół
vz
vvz
ad az
ad
ad ad
t [s]
Rys. 5 Przebiegi czasowe członu formującego wzorzec prędkości z niezróżnicowanymi
wartościami ograniczeniami przyspieszenia (ar = ah = ad )
Obwód twornika:
 bez uwzględniania nieliniowości przewodzenia przekształtnika
= - -
 z uwzględnieniem nieliniowości przewodzenia przekształtnika (wersja 1)
- - , e"

=
, < 0
 z uwzględnieniem nieliniowości przewodzenia przekształtnika (wersja 2)
- - , e" (" - - > 0

=
, < 0 (" - - d"
 9
Ogólna struktura programu symulacyjnego
Program symulacji
Dane : - czas symulacji,
- war. poczÄ…tkowy, opcje,
- uchwyt modelu.
Wywołanie solwera Ode
Model:
Ode: wywołania inicjacyjne
dxdt = model( t, x, p1, p2, p3= [] )
- Modelu:
dxdt = model( t, x, p1, p2, p3 )
Solver
- funkcji wyjścia:
Funkcja wyjścia wyfun:
wyj = wyfun( t, x init , p1, p2, p3 )
dane = []
Model:
danet = model( t, x, p1, p2, p3 =2 )
danet = model( t, x, p1, p2, p3 =2 )
dane = [dane; danet]
Nowa sekwencja
Ode: sekwencja wywołań z
danymi punktów obliczeniowych
Model:
[t,x]:
dxdt =model( t, x, p1, p2, p3= [] )
dxdt = model( t, x, p1, p2, p3 =[] )
Wyznaczenie i ocena punktu
rozw. numerycznego [t,x]
Ocena pozytywna ?
N
T
Ode: archiwizacja punktu [t,x] Funkcja wyjścia wyfun:
Model:
danet = model(t ,x, p1, p2, p3 =1 )
Wywołanie z danymi punktu [t,x]
danet = model( t, x, p1, p2, p3=1)
dane = [dane; danet]
wyj = wyfun( t, x, [] , p1, p2, p3=1 )
N
Koniec obliczeń ?
Funkcja
T
wyjścia
Ode: Wyjście  rozw. numeryczne [T, X]
Model
Odczyt archiwum funkcji wyjścia
Funkcja wyjścia wyfun:
dane = wyfun( [], [], 'danewy', [], [], [] )
Wyjście: dane
Odczyt archiwów lokalnych
Archiwizatory lokalne:
dane_loc = wyfun_loc( [], [], & , 'danewy' )
Wyjście: dane_loc
Wykresy, wydruki
Rys. 1.1-3 Schemat blokowy programu symulacyjnego
(kursywą oznaczono przykładowe instrukcje Matlaba)
 10