WFW 87.014
OEFENOPGAVEN M E N MET 'PC-MATLAB'
Verslag in het kader van een stage-opdracht. Sectie Regeltechniek Vakgroep WFW Faculteit W ~ r k ~ u i g ~ o u w ~ ~ l ~ ~ e Technische Universiteit Eindhoven
W.M.H. Vaassen Februari 1987
Prof.Dr.Ir. J.J. Kok dank ik voor het begeleiden van de stage.
- 1 -
. INHOUDSOPGAVE :
pag *
1
Inleiding
2
2
Wat is PC-Matlab?
3
3
Alvorens te starten met het PC-gebruik. 'I PC-DOS 2 Printers 3 Editors 4 'Onoplosbare' prohlenen
4
Leren werken met PC-Matlab. 1
2 3 4 5
6
De handleiding Diskettes Opstarten On-line Help ' .m-files' ' Diary-files'
Appendix : Oefenopdr~~h~en A1
Oefening 1: Toestandsbeschrijving van een massa-veersysteem.
A2
Oefening 2: Berekening van de exponent van een vierkante matrix.
A3
Oefening 3: Toestandsbeschrijving en Laplacetransformatie.
A4
Oefening 4: Optimale regeling en optimale toestansschatting van een lineair tijdskonstant systeem.
- 2 -
Hoofdstuk1
Inleidinq In dit stageverslag wordt een aantal mogelijkheden beschreven om met behulp van een Personal Computer (PC) praktische oefeningen uit te voeren. Het gaat hier speciaal over de toepassingsmogelijkheden van het programmapakket PC-Matlab in het kader van de kollegevakken Dinamika, Regelen, Fourlertechniek en stochastische signaalverwerking en Random Vibrations. De nadruk ligt in eerste instantie op de eerste twee genoemde vakken Het doel van het verslag is het weergeven van een globaal overzicht van onderwerpen die aan de orde komen bij het exerceren met PC-Matlab. Denk bijvoorbeeld aan het bestuderen van diverse handleidingen. In de appendix wordt een viertal oefeningen weergegeven die eventueel in de vorm van een opdracht en zonder beseleidinq gemaakt kunnen worden.
- 3 -
Hoofdstuk 2
Wat is PC-Matlab? 'Natlab' is een afkorting van 'matrix laboratory'. PC-Matlab is een voor PC's geschikt programmapakket. Het is afgeleid uit de projekten LINPACK en EISPACK en geschreven in de taal ' C ' . PC-Matlab is in hoofdzaak gericht op het uitvoeren van numerieke berekeningen met matrices. Deze matrices hoeven niet te worden gedimensioneerd! PC-Matlab werkt interaktief. Het gebruik ervan is relatief eenvoudig; men het programma opstart kan men matrices, dus ook skaìars, definieren. De resultaten van bewerkingen daarop zijn direkt zichtbaar. Het is tevens mogelijk routines aan te roepen en plotjes van bv. signalen te maken. Als
Men kan zelf routines schrijven en deze toevoegen aan de standaardroutines. Bovendien is er een 'HELP'-faciliteitdie on-line informatie geeft over standaard- gg zelfgeschreven routines! Echter, de werking van sommige programmaroutines is vrij traag en de input/output faciliteiten zijn beperkt.
- 4 -
Hoofdstuk 3
Alvorens te starten met het PC-aebruik.
I 3.1
Personal Computer Disk Operatina System.
Om met PC-Matlab te kunnen exerceren is het niet noodzakelijk een uitgebreide handleiding van het zgn. Disk Operating Systeem (PC-DOS) diepgaand t e bestuderen. Een overzicht van de meest elementaire handelingen is voldoende. Een voorbeeld van zo'n overzicht is de RC-Informatie AG-89 (25 blz.). Deze heeft tot doel een beginnende PC-gebruiker op weg te helpen. Tevens wordt enige achtergrondinformatie met betrekking tot PC's gegeven.
3.2
Printers
Verder is het nuttig te weten hoe men een bepaald type (merk) printer kan aansturen ten behoeve van het printen van een textfile van een programma o f het plotten van een plaatje. Het is aan te bevelen hier een korte handleiding voor te naken, temeer daar deze praktische problemen veel tijd kunnen kosten als men geen begeleiding heeft. In de handleiding hoeven slechts enkele veel gebruikte kombinaties van PC en printer aan de orde komen.
- 5 -
§ 3.3
De 'Editor'
Voor het schrijven van programmafiles is de beschikking over een eenvoudige 'editor' (wordprocessor, textverwerker) noodzakelijk. De DOS-line-editor EDLIN is te beperkt voor het schrijven van PC-Matlab programmafiles. De 'Turbo-Pascal' editor is uitgebreid genoeg en wegens zijn beperkte geheugengebruik geschikt om toe te passen in kombinatie met PC-Matlab. Een overzicht ( 1 pagina A 4 ) van de edit-kommando's is noodzakelijk en voldoende om 'Turbo' te kunnen gebruiken.
Q 3.4
'Onoplosbare' Problemen
Een beginnende PC-gebruiker heeft nog al eens te kampen met problemen waarvan de oorzaak hem geheel onbekend is. Voorbeeld: de komputer reageert nergens meer op; alleen uit- en weer inschakelen helpt nog! Dit kan het gevolg zijn van het niet op elkaar afgestemd zijn van apparatuur enfof software. Het is erg onbevredigend als er regelmatig onverwachte verschijnselen optreden terwijl de gebruiker niet weet of deze het gevolg zijn van systeemfouten ofwel van een fout die hij zelf, onbewust, steeds weer opnieuw maakt. Het is daarom gewenst dat hij de oorzaak van dit soort problemen kan achterhalen opdat hij niet steeds met dezelfde moeilijkheden gekonfronteerd wordt; in deze gevallen moet hij de hulp kunnen inroepen van iemand met een gedegen kennis op dit gebied.
- 6 -
Hoofdstuk 4
Leren werken met PC-Matlab.
3 4.1
De handleidinq
Bet is ~ a n z e ~ ~ ~ p r edat ~ eaen n d eerst de handleiding van PC-Ratlab bestudeert als men met het programma gaat werken. Dit houdt echter niet in dat de handleiàing in zijn geheel bestudeerd moet zijn voor dat aen kan starten! in eerste instantie is het voldoende als de opzet van de ~ ~ n d ~ e bekend i ~ i ~is; g
Hen kan de volgende indeling aanbrengen:
- Opstarten - A l g e ~ efunkties ~~
- Geavanceerde routines Het 2" en 3e deel bestaan ieder uit een 'Tutorial' en een 'Reference'. In de 'Tutorial' worden wat algenene toepassingsvoorbeelden behandeld terwijl ia de 'Reference' alle items op min o f neer alfabetische volgorde aan de orde komen. De hand~eidingheeft een omvang van ca. 340 pagina's formaat A5. Net is niet erg zinvol de hele handleiding te bestuderen en pas daarna met PC-Matlab te gaan werken. In tegendeel; door eerst te proberen en daarna eventueel de handleiding te raadplegen leert men veel sneller.
f 4.2
Diskettes
PC-Matlab kan wor6en gebruikt in de vorm van 2 diskettes. Op hen diskette de editor en een opstartprogramma. D i t is de staan o.a. het hoo~dprogranma~ zgn. 'Program-diskette'. Op de andere, de 'Utility-diskette',staan alle routines van PC-Matlab. Z e zijn gegroepeerd onder de sub-directory \MATLAB. In de root-directory staan de zelfgeschreven routines. Als er twee diskdrives beschikbaiir zijn, dan hoort de 'Program-diskette' in de A-drive en de 'Utility-diskette'in de B-drive.
- 7 -
$ 4.3
013starten
Het opstarten van PC-Natlab staat beschreven in het eerste hoofdstuk van de handleiding. Echter, de beschrijving is nog al onoverzichtelijk en bevat bovendien fouten. Hen kan beter een vervangende instruktie gebruiken: Om het programma te starten dient men een aantal kommando's te geven. Deze kommando's zijn voor het gemak geprogrammeerd in een 'batchfile' met de naam AUTOEXEC.BAT. De kommando's worden automatisch uitgevoerd als het Disk Operating Systeem wordt opgestart.
De inhoud van AUTOEXEC.BAT: timer/ s mgkeybrd int60h int61h path=a:\;b:\ set WTLABPATH=b:\matlab;b:\ b: echo off echo Start PC-MATLAB met het kommando: ' pcm ' echo of de editor met het kommando: * ed '
Na het draaien van AUTOEXEC.BAT is de default-drive B. Het heeft immers weinig nut om de default-drive A te laten zijn. Er worden daarop toch geen routines bewaard. Bovendien is er niet veel ruimte meer vrij op de 'Program-diskette'. Als men bv. een zgn. 'Diary-file' laat aanmaken, dan komt deze op de 'Utility-diskette'. Ook als men de editor start worden de geschreven texten automatisch op de 'Utility-diskette' bewaard.
- 8 -
$ 4.4
On-line Helu
In het hoofdprogramma van PC-Matlab is er 'on-line' informatie ter beschikking. 1
Er is een programma 'demo' dat toepassingen demonstreert van allerlei funkties .
2 Er kan een kommando 'help' gegeven worden. Een kombinatie 'help' 'item' is ook mogelijk. Zelfs 'help' 'help' geeft nuttige informatie.
Het is raadzaam eerst gebruik te maken van de help-faciliteit en pas in tweede instantie de handleiding te raadplegen.
4.5
'.m-files'
Hen kan zelf programma's schrijven in 2 vormen: 1
Function-file Deze bevat een aantal opdrachten en kan worden aangeroepen met 1 of meer argumenten. Variabelenamen die in zo'n programma gebruikt worden zijn lokaal, d.w.z. alleen binnen de routine van betekenis.
2
Command-file Ook dit type fife bevat een serie opdrachten, maar het laten werken van zo'n file heeft hetzelfde resultaat als het 'met de hand' uitvoeren van de opdrachten. Gebruikte variabelenamen hebben een globale betekenis.
Deze programmafiles heten '.m-files' omdat PC-Matlab alleen files met de extentie ''m' achter de naam als programma herkent.
- 9 -
3 4.6
'Diary-files'
Tijdens het editen van een .m-file is men niet meer binnen het hoofdprogramma van PC-Matlab. De on-line helpfaciliteit is dan niet meer ter beschikking! Om een statement uit te proberen moet men achtereenvolgens:
- het reeds geprogrammeerde deel opslaan - uit de editor gaan
-
het hoofdprogramma opstarten
- de gewenste tests uitvoeren
-
uit het hoofdprogramma gaan - de editor starten - de 'workfile'-naaminvoeren - het programmeren voortzetten Dit is niet alleen vervelend; de kans is groot dat men de testresultaten weer vergeten is voordat er weer geprogrammeerd kan worden! Bovendien kost het vaker verrichten van deze handelingen veel tijd. Het genoemde probleem is gedeeltelijk te verhelpen door bij het opstellen van een .m-file gebruik te maken van een 'diary-file'. Hierbij kan men de volgende werkwijze te hanteren:
-
-
Geef in PC-Matlab het kommando 'diary "naam" ' . Test alle kommando's die geprogrammeerd moeten worden. Geef het kommando 'diary off'. Verlaat PC-Matlab. Print de file "naam". Delete de file "naam" i.v.m. disketteruimte. Start de editor en schrijf de programmafile.
Als voorbeeld zijn vier oefenopdrachten opgesteld. Hierin wordt telkens eerst een inleiding gegeven waarna eventueel enkele handberekeningen worden gevraagd. Tot slot volgt een aantal opgaven die met PC-Matlab uitgevoerd moeten worden. Bij laatstgenoemde opgaven worden op sommige plaatsen aanwijzingen gegeven betreffende de te gebruiken Matlab-statements. Deze aanwijzingen zijn vet gedrukt. Ze worden voorafgegaan door een *. Ook de uitwerking is aan iedere opgave toegevoegd in de vorm van I f eventuele 'handberekeningen'
2 ) de letterlijke tekst van een .m-file 3)
numerieke en grafische uitvoer
Hierbij dient opgemerkt te worden dat de numerieke resultaten door PCMatlab in een minder kompakte vorm op het beeldscherm weergegeven worden. In dit verslag zijn ze echter qua vorm gemodificeerd ten behoeve van de leesbaarheid.
APPENDIX
A 2.5
Resultaten met
e1gdiag.m :
k o n d i t i e g @ t a l = 5.0000E+015
€2
v =
=
2.7183
o
O O O O
2.7183 O O O
1.0000
I . O000 -0.0000
O O
O
O O
1.0000 O
o O
o
U
O =
1 O O
o O
>
O 0 2.7183 O O
O 1 O 0 O
O U 1 O O
V(2,2)
anc = - 4 . 0 0 0 0 E - 0 1 6
O
O O O
1 O
O
o
O O
1
O
O
Q 2.7183 O
O O
o
o
O O O 1.0000 O
2.7183
O O O O 1.0000
A 3.2
Transformeer de bewegingsvergelijking tat een toestansbeschrijving ( z i e Regelen 11).
Kies als toestansvektor: &(t) =
Uitgangsvektor:
Ingangsvektor:
Y(t) -
[ i: ]
= q(t)
A 4.2
PC-MATLAB : Definieer de konstanten a en b. Opdracht 1 Stel de matrices A , B en C op. Bereken de regelbaarheidsmatrix en de rekonstrueerbaarheidsmatrix en bepaal of het systeem regelbaar enjof rekonstrueerbaar is.
*
ctrb
obsv
rank
Opdracht 2 Stel een tijdarray Ikofom) t 301 elementen.
op
dat gaat van O t / m 900 sekonden met
Genereer een ruisarray w m.b.v. t. Kies een normale kansdichtheidsverdeling met een standaarddeviatie van O .O5 m sek-2. Bereken m.b.v. de autokorrelatie van w de korrelatiematrix Vww van de vektor g(t). Plot w als funktie van t.
* *
t = [ 0:3ûû ] ' * 3 ; rand xlabel pause echo on
cov
echo o f f
plot
title
Opdracht 3 Bereken de eigenwaarden van A . Bereken en plot de responsie yft) op de ruis w(t) als u(t)
*
eig
5
O.
A
1.1
OEFENING
1
Toestandsbeschrijvinq van een massa-veersysteem.
11
Stel. de bewegingsvergelijking in y(t) van onderstaand massaveerxysteea met 1 graad van vrijheid op, uitgedrukt in wo en b Kies hierbij wo = J(k/ml en 5 = 2 J (m. k 1
1,
figuur 1.1
2)
Voer 3 variabelen in; xl(t) = y(t) x2(t) = Y W u(t) = F(t)
Voeg xl(t) en x2(t) samen in de toectandsvektor
x(t) =
[ :i:::J
A 1.2
Transformeer de bewegingsvergelijking m.b.v. de ingevoerde variabelen tot een svsteemverqeliikinq en een uitaanssveraeliikinq in de volgende vorm:
Stel de matrices A, b en c op.
A 1.3
3)
Stel i:
>
(Bovenkritische demping);
1
Bereken de eigenwaarden van de systeemmatrix A uitgedrukt in wo en S. Bepaal tevens de eigenvektoren van A uitgedrukt in de bijbehorende eigenwaarden.
4)
Stel O
<
i:
<
1
(Onderkritische demping);
Bereken de eigenwaarden van de systeemmatrix A . Bepaal tevens de eigenvektoren van A . Merk op dat de eigenwaarden en eigenvektoren komplex zijn!
5)
Stel i: = 1
(Kritische demping);
Bereken de eigenwaarden van de systeemmatrix A . Bepaal tevens de eigenvektoren van A . Merk op dat de eigenwaarden gelijk zijn en dat ook de eigenvektoren gelijk zijn; de matrix van eigenvektoren is niet van volle rang!
A 1.4
Voer de volaende ordrachten uit met PC-EVLTLAB:
Stel een tijdarray t op dat in 200 stappen het tijdsinterwal [ O , 6 ] sek. diskretiseert en een inputsignaal u(t) dat identiekgelijk aan O is.
*
t = (0:200)'*6/200;
u = zeros(t1;
Stel de systeemmatrix Al op. Bereken van A , de matrix V1 van de eigenvektoren en de diagonaalmatrix Dl van de eigenwaarden.
*
[ V1, D1
f = eig(A1)
Bereken de impulsresponsie van het systeem.
* *
i1 = impulse(
A l , b, c, 0, 1, t
1;
plot( t, il), title('1ipulsresponsie (Bovenkritisch gedempt)')
Bereken de stapresponsie van het systeem.
*
s l = step( A l , b, e, 0, 1, t
1;
A 1.5
We noemen x ( t ) een 'mode' van het systeem als i(t) = pex(t) en x(t) is reëel en # Q voor t>O. Een mode is de vrije responsie van het systeem voor t>O met begintoestand &(O) =
q , = aex waarbij y een reële eigenvektor van de systeemmatrix en
u
een reëel getal is.
Bepaal alle modes van het systeem. Bereken daartoe de vrije responsies van liet systeem met als begintoestand steeds een veelvoud van een reële eigenvektor. Kies dit veelvoud z.d.d. de beginverplaatsing 1 is.
Geval 2:
wo = 2
1
en 2, = -
;
2
Stel de systeemmatrix A2 op. Bereken van A2 de matrix V2 van de eigenvektoren en de diagonaalmatrix D2 van de eigenwaarden. Bereken de impulsresponsie van het systeem. Bereken de stapresponsie van het systeem-. Bepaal indien mogelijk de modes van het systeem.
A 1.6
Geval 3:
wo = 2
en
<
= 1 ;
Stel de systeemmatrix A3 op. Bereken van A3 de matrix V3 van de eigenvektoren en de diagonaalmatrix D3 van de eigenwaarden. Bereken de impulsresponsie van het systeem. Bereken de stapresponsie van het systeem. Hoeveel modes heeft dit systeem? Bepaal alle modes.
Uitwerking:
2)
Toestandsbeschrijving in expliciete vorm:
b = p ]
c = [ 1 ;
o1
a 1.8
3)
Eigenwaarden als Z
>
1:
Al
-
{
h2 = i
-z -
J ( 52
-
1 1
-5
J ( 22
-
1
1
Eigenvektoren:
4)
Eigenwaarden als
XI I ] en
+
[;I
z < 1:
Eigenvektoren:
5)
Eigenwaarden als 5 = 1:
Beide eigenvektoren zijn gelijk aan:
1
}.Wo
)'W*
A 1.9
PC-Matlab Commandfile:
echo on , cla
% MODES % % Commandfile als voorbeeld van d e uitwerking van oefening 1. 1. Laat dit programma draaien met '
1. 2 -De oefening betreft d e toestandsbeschrijving van een eenvoudig 1 massa-veersysteem met î graad v a n vrijheid. 1. -De eigenwaarden en elgenvektoren v a n het systeem worden berekend. 1. -Impuls- en stapresponsies en eigenbewegingen (modes) komen aan d e orde.
t = (0:200)'*6 /2O8
b = C O , 1 3 ' ,
:
u = zerosít):
c = C l , û I
1. Geval 1 tbovenkritisch gedempt)
C V I , D l 1 = eig(A1)
it = impulse( A l , b , c , O . 1 , t); echo o f f plot(t,ili title('1mpulsresponsie [bovenkritisch)') pause, echo on
A 9.10
sl = step( A l , b, c , O , 1 , t 1 ; echo o f f plot(t,sl) title('Stapresp0nsie (bovenkritisch gedempt)') pause, echo on
xo11 = Vl(:,l) , xo12 = V 1 ( : , 2 ) / V 1 ( 1 , 2 1 rll = lsim( A l , b, c , O , u , t , xLìlll; r12 = Isimí A l , b. c, O , u, t , ~ 0 1 2 ) ; echo o f f ploti t , C rll, r12 1 title('Modes (bovenkritisch gedempt)') pause, echo on
1
2 Geval 2 (onderkritisch gedempt) A 2 = C O , 1; - 4 , - 2
1
1 V 2 , O 2 1 = eig(A2)
i2 = impulsel A Z , b, c, O, 1 , t);
echo o f f plot(t,i2) titlet'Impulsresponsie (onderkritisch)') pause, echo on s2
= step( A 2 , b, c, O , 1 , t
1;
echo o f f plot(t,s2f titleí'Stapresponsie (onderkritisch gedempt)') pause, echo on
2 Geval 3 (kritisch gedempt)
A3 = c] O , 1 ; - 4 ,
-4
E V 3 , D3 3 = eig(A3)
I
A 1.11
i3 = impulse( A 3 , b , c , O , 1 , t);
echo o f f plotít,i3) title('Impu1sresponsie (kritisch gedempt)') pause, echo on s3 = step( A 3 , b , c, O , 1 , t 1 : echo o f f plot(t,s3) title('Stapresp0nsie (kritisch gedempt)') pause, echo on x03 = V 3 ( : , 1 ) + V 3 ( : , 2 ) ; ~ 0 3= x 0 3 / ~ 0 3 ( 1 1 ra = lsPm( A 3 , b , C , O . U . t , X Q 3 1 ;
echo o f f plot( t, r 3 title('1 Mode (kritisch gedempt)') pause, echo on
1 samengevat
echo o f f plot ( t, C il, i2, i3 1 1 title 'Impulsresponsies') pause plot ( t , E S I . s2, s 3 1 1 title 'Stapresponsies') pause plot ( t , E rll, r12, r 3 I i title('2 Bovenkritisch- en 1 kritisch gedempte mode') pause, echo on who 2 EINDE PROGRAMMA 'MODES'
A 1.12
Resultaten:
b =
0
c =
1
1
0
X Geval I (bovenkritisch gedempt)
AI =
Vl =
1
0 -4
-8
1.0000
-0.1340 1 * 0000
-0.5355
DI
o
-0.5359
-7.4641
O
,>:
.?:
J
_ .>' : -!
::
'.i
::i:
A 1.13
xoli =
2.0000 -0.5359
xo12
=
1.0000 -7.4641
A 1.14
X G e v a l 2 fonderkritisch gedempt)
A2 =
V2 =
O
1
-4
-2
-0.2500
-
0.433Oi
1.0000
O2 =
-1.
O000 + 1.7321i
o
- 0 . 2 5 0 0 + 0.4330i 1 .O000
0
-1.0000
-
1.7321i
1
A 1.15
1 Gewal A3 =
3
(kritisch gedempt)
0
1
-4
-4
V3
=
-0.5000
D3
=
-2
a
O
-2
-0.5000 1 .o000
1 .oooo
A 1.16
x03
=
1.0000 -2.0000
Z samengevat
A 1.17
four variables are:
eps, pi, inf, nan,
...
Al A2
D3 VI
C
r12
i1
i-3
U
A3
o1
v2
i2
SI
v3
i3
s2
D2
b
rt 1
s3
xO11 XOI 2 x03
t
A 2.1
OEFENING 2
Berekenincr van de exponent van een vierkante matrix.
Het berekenen van een overgangsmatrix van een tijdskonstant systeem met systeemmatrix A kan op verschillende manieren gebeuren (Regelen I I ) : machtreeksbenadesing, diagonalisatie, inverse Laplace-trans£ormatie. In deze oefening komen de eerste 2 methoden aan de orde. Beschouw daartoe de berekening van eA waarbij A een vierkante matrix is.
Machtreeks:
-
Kopieer het programma \matlab\expm2.m naar de hoofddirectory b:\ onder de naam mreeks-m Doe dit onder PC-DOS.
- Ga
na wat de funktie van het PC-Watlab kommando 'norm(A,1)' is.
- Wijzig mreeks.81
;
Het programma moet een array opsteilen waarin het ie element de norm van het ie inkrement van de (geschaalde) reeks voorstelt. De andere funkties moeten behouden blijven. De aanroep moet als volgt zijn: [€,inorail = mreeks(A) .
*
inorm = [ inorm, norm(F) 3
A
2.2
Diauonalisatie:
-
Kopieer het programma \matlab\expm3.i naar de hoofddirectory b:\ onder de naam eigdiag-n. Doe dit onder PC-DOS.
- Wijzig eigdiag.i
;
Het programma moet ook het konditiegetal van de matrix van eigenvektoren uitvoeren naar het beeldscherm.
PC-Matlab:
Bepaal met mreeks en met eigdiag de exponent van diverse vierkante matrices: algemene symmetrische met gelijke eigenwaarden niet-symmetrische met gelijke eigenwaarden
A 2.3
PC-flatlab Function-files mreeks.m en eigdia9.m :
t
function CE,inorml = mreekst.4)
Z E = Matrix exponential of A via Taylor series. Z inorm = C norm(fl),norm/F2), 3 IFi = scaled increment!
...
2 Synopsis
íE,inorml = mreeks(A1
Z Scale A by power of 2 so that its norm is < 1 / 2 s = round(log(norm(A,?)1/logt2)+1.5); if s < O , s = O ; end P. = b J 2 - s ;
.
Z Taylor series for expiAl E = O*A; F = eye(A1; k = I ;
inorm = normtfl; while norm(E+F-E,î) > O inorm = Einorm,norm(Ffl; E = E + F; F = A*F/k; k = k+l; end
1 Undo scaling by repeated squaring for k = 1 : s ' E = €*E; end
2
function E = eigdiaglA)
Z Matrix exponential via eigenvectors and eigenvalues
1 Print het konditiegetal van d e matrix van eigenvektoren Z Synopsis
E = eigdiagíA)
tv,dl = eig(A1; konditiegetal = cond(v) E = v * diag(exp(diagldl)) / v;
A 2.4
Resultaten ntet mreeks.m :
A =
I O O O O
>
[E?,inorml=mreeksfA)
EI =
2.7183 O
2.7183 2.7183
O O
O O O
o
inorm =
O O 2.7183
O O O
U
2.7183
O O O O
O
O
2.7183
O. 2023
0.0189
0.0011
o . O000
0.0000
o . O000
o . O000
o . OOOO
-0.6941
-1.7244
-2.9685
-4.3656
-12.7281
-14.5864
-16.4902
Columns 1 through 7 1 .o000
1.0000
o. O000
Columns 8 through 13
o. O000
o
*
0000
> log10íinorm) ans =
Columns
O
1
through 7
O
Columns 8 through 13 -7.4861
-9.1704
-10.9205
-5.8792
A 3.1
OEFENING 3
Toestandsbeschriivinff en LaPlacetransformatie.
De beweginqsvergelijking in q(t) van onderstaand massaveersysteem luidt:
figuur 3. I
Bepaal met behulp van de methode van Lagrange de matrices M, B en K in de bewegingsvergelijking als funktie van de parameters m, b, en k.
A 3.3
PC-Hatlab:
Stel de matrices M, B en K op voor bepaalde zelfgekozen parameters.
Bereken de matrices A , 8 , C en D van de toestandsbeschrijving uit de matrices Pi, B en K van de bewegingsvergelijking
Bereken de eigenwaarden van A.
*
eig
Bereken de overdrachtsfunkties
*
-Y ( s )
en
-Y ( s ) .
ss2tf
Bereken de polen van de overdrachtsfunkties en vergelijk deze met de eigenwaarden van A.
*
roots
Maak Bode-plots van de overdrachtsfunkties t.a.v. zowel de amplitude als de fase.
*
logspace
bode
loglog
title
semilog
Maak Nyquist-plots van de overdrachtsfunkties.
*
axis('sguare')
polar
axis( ' norinal ' 1
A 3.4
P C - M a t l a b Command-file mvsyst.m :
echo o f f , clear, cla
i! 'nvsyst.m' stelt de matrices M , E , en K op voor een massa-veersysteem
stelt d e matrices A , E , C en R op voor een toestandsbeschrijving berekent d e eigenwaarden van de systeemmatrix berekent d e overdrachtsfunkties maakt Bode- Eon Nyquist-plots van d e overdrachtsfunkties
2 1 2 i! i!
1 Re invoer wordt opgevraagd. m l = input('rn.i = ' I ; m2 = inputí'rn2 'I: bl =
inputí'bl = ' 1 ; b2 = inputl'b2 'I; k3 = inprnt('b3 = ' I ; kl = inputf'k? = ' 1 ; k 2 = input('k2 = ' 1 ; k 3 = input('k3 = ' 1 ;
M = lm1,0;~,m21 K = [kl+k3,-k3;-k3,k2+k33 B = lbl+b3,-b3;-k3,b2ib31 pause
A B C D
= = = =
Izeros(2),eye Czeros(2i;inv Ceyeí21,zeros CzesQsf2)f;
pause
A 3.5
CNum1,denll = s s Z t f í A , B , C , D , l ) tNurnZ,den21 = ssPtf(A,B,C,D,Z) pause
eigenwaarden = eigíA) polen = rootsfden11 pause
opnieuw = 1 ; while opnieuw > 0 , 10 = input('6eef de logaritme van d e laagste hoekfrekwentie ' 1 : up = input('Feef d e logaritme van d e hoogste hoekfrekwentie ' 1 ; N = inputí'Aantal plotpunten * I : w=logspaceílo,up,N);
iog90g(w,MagI)
titlel'Magnitude v a n overdracht v a n F 1 ' ) pause semilogxíw,Phaseí) titlef'Fase van overdracht van F l ' ) pause, grid, pause loglog(w,Mag2)
title('Magnitude van overdracht van F Z ' ) pause semilogxíw,fhase2) title('face van overdracht van F 2 ' ) pause, grid, pause
A 3.6
axis('square') polaríPhase1,Magl) titleí'ûverdracht van F l ' ) pause, g r i d , pause poi.ar(Phase2,MagZl title('0vesdracht van F 2 ' ) pause, g r i d , pause axis('normal'1 opnieuw = input('frekwentiedomein aanpassen ( 1 1 , anders 1 0 ) ' 1 : end
cla echo on who
x
einde programma
A 3.7
Resultaten:
ml = 3 m2 = 1 bl = 0 . 3 b2 = 0.8 b3 0.2 kl = 9 k2 = 5
k3 = 2
M =
3
O
O
I i
K =
11
-2 7
-2
B =
0.5000 -0.2000
A =
n
I
O
O
1 O000
O
O
O O . 6667 -7.0800
O -0 1667 o . 2000
1.0000 O . 0667 - 1 * 0000
-3.6657
2.0000
O O . 3333
O O O
O
1 . o000
8 =
c =
-0 2000 1 O000
0
1 O
O 1
O O
s
e
O O
i
A 3.8
o. O 0 0 0 o. o 0 0 0
O. 3333
O
o. O000
O. 3333 O . 0667
2.3333 O. 6667
den1 =
1.0000
1.1667
10.8200
4.5667
24.3333
Num2 =
O O
o . O000 o . O000
-0.0000 I . O000
O . 0667
O . 6667
O. 1667
3.6667
2.0000
1.1667
1 O. 8200
4.5567
24.3333
Numí =
den2 =
O
eigenwaarden =
polen =
-0.4969 -0.4969 -0.0864 -0.0864
+ 2.6609i
- 2.66091 + 1.82031 - 1.82031
- 0 . 4 9 6 3 + 2.6603i - 0 . 4 9 6 9 - 2.6609i -0.0864 + 1.8203i -0.0864 - 1.82031
A 3.9
..
. <
.. ..
... .
.,
.. .
... ..
..., ..,, ..,. .... .... . . . ..
.
..
.. .. .. .. ,.
I
1
A 3.10
.
_ "
.
.
.
A 3.11
Your v a r i a b l e s are:
eps, p i , i n f , nan,
...
A B
den1
m2
c
N Numl Num2
den2 elgenwaarden
opnieuw polen
D
Phase1
kl
fc
Phase2
UP
M
W
bl
k2
Mag1 Mag2
b2 b3
k3 lo ml
A 4.1
Oefenincr 4
Optimale reselinu en optimale toestandsschattins van een lineair tiidskonstant systeem. Beschouw het voorbeeld van de bathyscaaf in Regelen I1 pag. 3.122. Het gedrag van de bathyscaaf kan weergegeven worden met een systeemmodel als in figuur 4.1 .
figuur 4.1
O ]
C=[1
D = [ O ]
h(t) is de vertikale verplaatsing t.o.v. de duikdiepte 1000 m.
w(t) is de systeemruis
(mi
sek-')
.
u(t) is de verandering van de massa van de bathyscaaf t.o.v. de massa waarbij deze in evenwicht is op 1000 m. diepte. y(t) is de enige meetbare toestandsgrootheid h(t). a2 = l.01elû'4
seke2
m kg-1sek-2 b = ~3.4175-10-~
A
4.3
Veronderstel dat het systeem gestabiliseerd kan worden door terugkoppeling van x ( t ) . Dit kan weergegeven worden door het model in figuur 4.2 .
I
figuur 4.2
Er geldt nu:
Opdracht 4 Bepaal L, 2.d.d. de eigenwaarden van (A-BeL,) op -0.004 en -0.006 sek-l komen te liggen (zwakke terugkoppeling). Bepaal L, z.d.d. de eigenwaarden van (A-BeL,) op -1 .O04 en -1 .O06 sek-l komen te liggen (sterke terugkoppeling). Bereken g,(t) en gs(t) als responsie op w(t) bij terugkoppeling van de toestand met resp. L, en LS. Bereken ook de regelsignalen u,(t)
en us(t).
Plot samen de verplaatsingen als funktie van t bij zwakke resp. sterke terugkoppeling. Doe dit ook met de snelheden en de regelsignalen.
* * 2
*
1
Lz = place( A, B, [ -0.004, -0.006 xz = lsim( A-B*Lz, [ O 1 I ' , eye(21, zeros(2,1), wf t ) UZ = -B"Ez*xz' plot( t, xz(:tl), t, xs(:,l) 1
A 4.4
Merk op dat E[ u2(t) 1 groter is naarmate de terugkoppeling sterker is. Er bestaat een lineaire toestandsterugkoppeling Lo z . d . d . een gekozen minimaal is. kwadratisch kriterium J = E[ gT (t)eQeg(t) + u(t)*R*u(t) Hierin is Q semi positief definiet en R>O.
1
Opdracht 5 Stel een matrix (1 en een weegfaktor R op z.d.d. de snelheid van de bathyscaaf niet gewogen wordt terwijl 1 m. afwijking even zwaar meeweegt als een massaverandering van 1 kg. Bereken de optimale terugkoppelmatrix Lo en de regelpolen feigenwaardenvan (A-B*Lo)). Bereken &,(t) met Lo.
als responsie op w(t) bij terugkoppeling van de toestand
Bereken ook het regelsignaal u o ( t ) . Plot (een voor een) verplaatsing, snelheid en regelsignaal als funktie van t bij terugkoppeling met Lo. Bereken de standaarddeviaties van verplaatsing, snelheid en regelsignaal voor de drie gevallen zwakke- sterke- en optimale terugkoppeling.
*
lqr
SxzSxsSxo = std( [ XZ, x s f xo
3
1
A 4.5
Omdat de gehele toestandsvektor in werkelijkheid niet beschikbaar is zal een schatting gebruikt moeten worden. Deze schatting wordt gegenereerd met een waarnemer. De geschatte toestand noemen we &{t). De inputsignalen van de waarnemer zijn u(t) en y,(t). ym(t) is het met ruis gestoorde uitgangssignaal: y(t)+v(t). v(t) wordt aangenomen als een normaal verdeeld witte-ruissignaal.
..
..
figuur 4 . 3
Neem voorlopig aan dat v(t)
o O;
A 4.6
Opdracht 6 Bereken Ks z.d.d. de eigenwaarden van (A-Ks*C) op -0.500 en -1,000sek" komen te liggen. Bereken hs(t)als responsie op de input van de waarnemer. Kies als input yo(t) (komt overeen met k(t) in figuur 4 . 3 ) en uo(t). Plot samen de verplaatsing en de schatting daarvan. Doe dit ook met de snelheid.
Er wordt nu verondersteld dat er wel een meetruissignaal aanwezig is.
Opdracht 7 Stel een array v op dat het normaal verdeelde ruissignaal v(t) met standaarddeviatie 5 m. voorstelt en bereken de autokorrelatie Vw ( = ECv2(t)l 1. Bereken -Xwsv(t)als responsie op de input van de waarnemer. Kies als input yo(t) en uo(t). Plot samen de verplaatsing en de schatting daarvan. Doe dit ook met de snelheid.
Merk op dat de schattingen geen betekenis hebben door de veel te sterke terugkoppeling van de innovatie E(tf.
A 4.7
Opdracht 8 Bepaal de optimale terugkoppeling KO van (xtt)-xw(t)1 I. van EC fx(t)-s(t) Maak hierbij gebruik van Vww en Vw.
E(t)
t.a.v. het minimaliseren
Bereken de waarnemerpolen (eigenwaarden van (A-Ko*C)). Bereken &(t) als responsie op de input van de waarnemer. Kies als input y,(t) en u o ( t ) . Plot samen de verplaatsing en de schatting daarvan. Doe dit ook met de snelheid.
*
lqe
A 4.8
Als de geschatte toestand &(t) wordt teruggekoppeld i.p.v. de echte toestand x(t) dan zijn systeem en waarnemer gekoppeld. Dit kan worden weergegeven met het model in fig. 4.4.
I
figuur 4.4
De toestandsvergelijking van het gekoppelde systeem:
A 4.9
Opdracht 9 Stel de systeemmatrix A2 en de ingangsmatrix B2 van het gekoppelde systeem op. Beperk in dit geval de ingang tot
Bepaal de polen van het gekoppelde systeem (eigenwaarden van A Z ) . Vegelijk de polen met de regelpolen en de waarnemerpolen.
Qpdracht 10 Bereken x2(t) als responsie van het gekoppelde systeem OP w(t) en v(t) m.b.v. Al en B2. Vergelijk de verplaatsing en de schatting daarvan in een plotje. Doe dit ook met de snelheid. Vergelijk de verplaatsing bij terugkoppeling van de toestand met de verplaatsing bij terugkoppeling van de optimale schatting van de toestand. Doe dit ook met de snelheid.
A 4.10
PC-Matlab Commandfile:
echo on , clear , cla
1 2 1 2 1
BATHYSCA.M Optimale regeling van de duikdiepte van d e bathyscaaf (Regelen 1 1 ) . Berekening van een optimale statische toestands-terugkoppeling en een optimale waarnemer.
1 ENKELE KONSTANTEN a=cqrt(l.Ole-4),b=8.4i75e-4 pause cla
1 Opdracht 1 A = [ û , 1 ; a*a,O 1 , pause ela
B = [O;bl, C = [ 1 , 0 3
Controllability = ctrbtA,B), Controllabalityrang = rank(Controllabi1ity) pause cla Observability = obsv(A,C), Observabilityrang = rank(0bservability) pause cla
1 Opdracht 2 t=[0:3003'*3; rand('norma1') w = 0.05*rand(t) ; Vww = [ O , O ; O , cov(w) pause
1
1 Systeemruis m/sec2 2 Autocorrelatie van d e systeemruis
A 4.11
echo o f f plot ( t ,w) title('Systeemruis wit) pause, echo o n , cla
Im/secZJ')
Z Opdracht 3 eigenwaarden = eigtA) y = lsim(A, E O , 1 I ' , C, O , w , ti; pause echo o f f plotlt,y) titlei'Ruisrecponsie zonder regeling') pause, echo on, cla
Z Opdracht 4 Lt = place(A, B , E - 0 . 0 0 4 , -0.006 1 1 Z Zwakke toestandsterugkappeling L s = place(A, 6 , E - 0 . 1 0 4 , -0.106 1 X Sterke toestandsterugkoppeling xz = lsim(A-B*Lz, t O , 1 I ' , eyel21, zeresiZ,f), w, t); x s = lsim(A-B*Ls, [ O , 1 J', e y e ( Z f , zeros(Z,t), w, tl; u2 = -Lz*xz' ; us = -Ls*xc' ; pause
1
echo o f f plot( t, xz(:,l) , t , xs(:,ll) title('Verp1aatsing bij zwakke en sterke terugkoppeling') pause plot( t , xzl:,Z) , t , xs(:,21) title('Sne1heid bij zwakke en sterke terugkoppeling') pause plot( t ' UZ, t , us 1 title('Regelsignaa1 bij zwakke en sterke teTUgkoppeì.ing') pause, echo on, cla
A 4.12
1 Opdracht 5
a = [ 1 , O ; O, O ]
R =
2
1
l m
t.t
1kg
pause cla L o = l q r ( A , B , Q , R ) , regelpolen = eig(A-E*Lo) 2 Optimale toest.terugkoppeling
xo = lsim( A-B*Lo, C O, í J ' , eye(21, Z@rOS(P,f), w a t ) ; u0 = - L o * x o ' pause
;
echo o f f plotf t, x o f : , l I 1 title4'Verplaatsing b i j optimale terugkoppeling van d e toestand') pause plot( t, x o ( : , 2 ) title('Sne1heid bij optimale terugkoppeling van d e toestand') pauss plot( t , u0 1 title('Regelsignaa1 b i j optimale terugkoppeling van d e toestand') pause, e c h o on, cla
1
1 Standaarddeviaties van x en u b i j Zwakke, Sterke en Optimale terugkoppeling SxzSxsSxo = stdl C xz , SuzSusSuo = std( [ U Z ' , pause cla
XS
,
US',
XO
UO'
1 1 1
3
A 4.13
Z TOESTANDSSCHATTING m.b.v.
een WAARNEMER
1. Opdracht 6
Ks =
(
1 Sterke terugkoppeling
place( A',C', C -0.500 , -1.000 3 1 1 '
Bw = E Ks*C , B I ; xws = lsim( A-Ks*C, Bw, eye(21, meros(2,3ì, pause
XO
Uo'3, t);
I
echo o f f plot( t. xo(:,I) ' t , xws(:,1)) titïe('Verp1aatsing en schatting ZONDER meetruis') xPabeP('SYSTEEM en W A A R N E M E R WIET gekoppeld') pause plot( t, xo(:,2) ' t , xws(:,2)) title('Sne1heid en schatting ZONDER meetruis') xlabel('SYSTEEM en W A A R N E M E R NIET gekoppeld') pause, echo on, cla
2 Opdracht 7
z
v = S*rand(t); v v v = cov(vl
Z
Meetruis; standaarddeviatie 5 m op 1 km diepte Autocorrelatie van d e meetruis
Bw = I Ks*C , EL , K S I ; xwsv = lsimi A-Ks*C, Bw, eye(21, zeros42,4), pause echo o f $ plot( t , xo(:,l) , t , xwsv(:,.1)) title('Verp1aatsing en schatting MET M E E T R U I S ' ) xlabel('SYSTEEM en WAARNEMER NIET gekoppeld') pause plot( t, xo(:,2) t , xwsw(:,2)1 title('Sne1heid en schatting MET MEETRUIS') xlabel('SYSTEEM en WAARNEMER MIET gekoppeld') pause, echo on, cla 1
Xo
,
UO',
v 1 , t);
A 4.14
l Opdracht
6
K O = lqe( A , eye(21, C, Vww, V v v 1 waarnemerpolen = eig( A - Ko*C 1 BW = C K o * C , B , KO
l Optimaal filter
1;
xwo = lsim( A-Ko*C, Bw, eye(21, zeros(2,4ì , C xo , uo', v pause
I , t);
echo o f f plot( t , x0(:,1) , t , xwo(:,I)) titlel'Verplaatsing en OPTIMALE schatting MET MEETRUIS') xlabel('SYSTEEí4 en OPTIMALE WAARNEMER NIET gekoppeld') pause plot( t , X O ( : , 2 ) , t , X W Q ( : , 2 1 ] title('Sne1heid en OPTIMALE schatting MET MEETRUIS') xlabel('SYSTEEM en OPTIMALE WAARNEMER NIET gekoppeld') pause, echo o n , c l a
2 KOPPELEN VAN DE OPTIMALE WAARNEMER AAN HET TE REGELEN SYSTEEM
l Opdracht 9 Z
Matrices van het samengestelde systeem
AZ = [ A , -B*LO 82 = [ C O , 1 pause cla
1'
; KO*C,
A-B*LO-KO*C .I , zerosi2,I) : zerosi2,l). KO I
regelpolen,waarnemerpolen gekoppeldepolen = eigfA2I pause cla
A 4.15
Z Opdracht 10 x2 = lsim( A2, 62, eye(41, zeros(4,2), C w , v 1, t); echo o f f plot( t , x2(:,1) , title('Verp1aatsing xlabeli'SYSTEEM en pause plot( t , x2(:,2) , title('Sne1heid en xlabelf'SYSTEEM en pause
t , x2(:,3)) en schatting') OPTIMALE WAARNEMER gekoppeld') t , x2(:,4)) schatting') OPTIMALE WAARNEMER gekoppeld')
p l o t I t , XQ(:,!l
, t , X2[:,f)) title 'Verplaatsing b i j terugkoppeling van TOESTAND resp. SCHATTING') pause plot ( t , xo(:,2) ' t , x2(:,2)1 title 'Snelheid bij terugkoppeling van TOESTAND resp. SCHATTING') pause echo on, cla
Z EINDE PROGRAMMA 'bathyscaaf' who
A 4.16
Resultaten:
a = 0.0100 b = 8.4175E-004
1 Opdracht A =
1
O
t.
B =
1.OE-003
1
c =
uuoo O
0.0001
*
O 0.8418
O
Controllability = $.OE-003
*
O
0.8418
Controllabilityrang = 2
Observability = 4 O Observabilityrang = 2
X Opdracht 2 vww = o O
O O . 0027
O
1
8.8418 O
A 4.17
I
1 Opdracht 3
eigenwaarden =
0.0100 -0.0100
7. Opdracht 4
Lz =
0.1465
11.8600
Ls =
13.2165
269.4802
A 4.18
.
...
.
. i
A 4.19
1
Opdracht 5
Q =
I O
O O
R = l LO =
1.1272
regelpolen =
51.1507 -0.0218 + 0.01931 - 0 . 0 2 1 8 - 0,0193i
A 4.20
SxzSxsSxo =
144.7567
SUZSUSSUO =
24.0531
O . 6048
1.1778
36.6262
O . 1319
10 - 0 5 2 5
17.3530
O . 2552
A 4.21
1 Opdracht Ks =
6
1.5000 0.5001
A 4.22
S Opdracht 7
vvv = 24.2747
A 4.23
1
Opdracht 8
KO =
0.1460 0.0107
waarnemerpolen
=
-0.0730 + 0.0723i -0.0730 - 0.0723i
A 4.24
2
Opdracht 9
A2 =
O
1.0000 O O O
0.0001 0.1460 0.0107
O
82 =
-0.1460 -0.0115
0
-0.0436 1.0000 -0.0436
O O 0.1460 0.0107
1 .o000 O O
regelpolen
0
-0.0003
=
-0.0218 + 0.0193i -0.0218 - 0.0133i
waarnemerpolen =
-0.0730 + 0.0723i -0.0730 - 0.07231
gekoppeldepolen
= e i 9 (A21
gekoppeldepolen
=
X Opdrach t 10
-.,.:.-.
i .
j
ii.
-0.0730 + 0.0723i -0.0730 0.0723i -0.0218 + 0.0193i -0.8218 - 0.0193i
-
A 4.25
A 4.26
Your variables are: A A2 8 82 8W
C Controllability Controllabilityrang KO
eps, pi, inf, nan,
...
a b eigenwaarden gekoppeldepolen regelpolen t U0
us UZ
KS
V
Lo
W
Ls
waarnemerpolen
LZ
x2
Observability Observabilityrang
xs
Q
R
XO
xwo xws
xws v
suzsussuo sxZsxssx0
XZ
vvv
Y
VWW
