Oefenopgaven maken met 'PC-Matlab': verslag in het kader van een stage-opdracht

Citation for published version (APA):

Vaassen, W. M. H. (1987). Oefenopgaven maken met 'PC-Matlab': verslag in het kader van een stage-opdracht. (DCT rapporten; Vol. 1987.014). Technische Universiteit Eindhoven.

Document status and date: Gepubliceerd: 01/01/1987 Document Version:

Uitgevers PDF, ook bekend als Version of Record Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at:

openaccess@tue.nl

providing details and we will investigate your claim.

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

- 1 -

. INHOUDSOPGAVE :

1 Inleiding

2 Wat is PC-Matlab?

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 De handleiding 2 Diskettes 3 Opstarten 4 On-line Help 5 ' .m-files' 6 ' 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. pag * 2 3 A4 Oefening 4:

Optimale regeling en optimale toestansschatting van een

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;

A l s 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.

Men kan zelf routines schrijven en deze toevoegen aan de standaardroutines.

Bovendien is er een 'HELP'-faciliteit die 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.

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 A4) 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

Hoofdstuk 4

Leren werken met PC-Matlab.

3 4.1 De handleidinq

Bet is ~ a n z e ~ ~ ~ p r e ~ e n d dat aen 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 i ~ i ~ g bekend is;

Hen kan de volgende indeling aanbrengen:

-

Opstarten

- A l g e ~ e ~ ~ funkties

-

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~eiding heeft 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

staan o.a. het hoo~dprogranma~ de editor en een opstartprogramma. D i t is de

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

- 7 -

$ 4.3 013s tarten

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 'batch-

file' 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 int6 1 h 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

$ 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

f unkties

.

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'-naam invoeren

- 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.

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 PC-

Matlab in een minder kompakte vorm op het beeldscherm weergegeven

worden. In dit verslag zijn ze echter qua vorm gemodificeerd ten

Resultaten met e 1 g d i a g . m : k o n d i t i e g @ t a l = 5 . 0 0 0 0 E + 0 1 5 €2 = 2 . 7 1 8 3 o O O O O 2 . 7 1 8 3 0 o o O O 2 . 7 1 8 3 Q O O O O 2 . 7 1 8 3 O O O O O 2 . 7 1 8 3 v = 1 . 0 0 0 0 I . O000 O O O o - 0 . 0 0 0 0 O O O O O 1 . 0 0 0 0 O O O O O 1 . 0 0 0 0 O U o O O 1 . 0 0 0 0 O = 1 O O O O O 1 U O o O O 1 O O o 0 O 1 O O O O O

1

> V ( 2 , 2 ) anc = - 4 . 0 0 0 0 E - 0 1 6

A 3 . 2

Transformeer de bewegingsvergelijking tat een toestansbeschrijving ( z i e Regelen

11).

[

i:::

]

Kies als toestansvektor: &(t) =

Y(t) =

-

q(t)

-

Uitgangsvektor:

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 op dat gaat van O t / m 900 sekonden met

301 elementen.

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.

*

*

xlabel pause

t = [ 0:3ûû ] ' * 3 ; rand

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) 5 O.

A

1 . 1

OEFENING

1

Toestandsbeschrijvinq van een massa-veersysteem.

1 1

2)

Stel. de bewegingsvergelijking in y(t) van onderstaand

massaveerxysteea met 1 graad van vrijheid op, uitgedrukt in

wo en

1,

Kies hierbij

wo

= J(k/ml en 5 = b

2 J (m. k

1

figuur 1.1

Voer 3 variabelen in; xl(t) = y(t)

x2(t) =

Y W

u(t) = F(t)

[

:i:::

J

Transformeer de bewegingsvergelijking m.b.v. de ingevoerde

variabelen tot een svsteemverqeliikinq en een uitaanssveraeliikinq in de volgende vorm:

A 1 . 3

3 ) Stel i:

>

1 (Bovenkritische demping);

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!

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 identiek-

gelijk aan O is.

*

t = (0:200)'*6/200; u = zeros(t1;

Stel de systeemmatrix Al op.

Bereken van A , de matrix V 1 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.

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.

1

2

Geval 2: wo = 2 en 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-.

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.

Uitwerking:

2) Toestandsbeschrijving in expliciete vorm:

3 ) Eigenwaarden als Z

>

1 : A l

-

-

{

-z

-

J ( 52

-

1 1 } . W o h2 = i -5

+

J ( 22

-

1

1

)'W* Eigenvektoren: 4) Eigenwaarden als

z

<

1:

Eigenvektoren: 5 ) Eigenwaarden als 5 = 1: 1 XI I ] en

[;I

A 1 . 9

PC-Matlab Commandfile:

echo on , cla

% M O D E S

%

% Commandfile als voorbeeld van d e uitwerking van oefening 1.

1. Laat dit programma draaien met ' < C t r l > < P r t S c 9 ' a a n ,

Z mits er een printer aan de PC gekoppeld is.

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 ( m o d e s ) komen aan d e orde.

t = (0:200)'*6 /2O8 : u = zerosít):

b = C O , 1 3 ' , 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); e c h o o f f plot(t,ili title('1mpulsresponsie [bovenkritisch)') pause, echo on

sl = step( A l , b, c , O , 1 , t 1 ;

echo o f f

title('Stapresp0nsie (bovenkritisch gedempt)') pause, e c h o on p l o t ( t , s l ) xo11 = V l ( : , l ) , xo12 = V 1 ( : , 2 ) / V 1 ( 1 , 2 1 rll = lsim( A l , b , c , O , u , t , xLìlll; r 1 2 = Isimí A l , b. c, O , u, t , ~ 0 1 2 ) ; echo o f f ploti t , C r l l , r 1 2 1

1

title('Modes (bovenkritisch gedempt)') pause, e c h o on 2 G e v a l 2 (onderkritisch gedempt) A2 = C O , 1; - 4 , - 2 1 1 V 2 , O2 1 = eig(A2) i2 = impulsel A Z , b, c, O, 1 , t); echo o f f p l o t ( t , i 2 ) 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 I

A 1 . 1 1

i3 = impulse( A 3 , b , c , O , 1 , t); echo o f f

p l o t í t , i 3 )

title('Impu1sresponsie (kritisch gedempt)') pause, e c h o on

s3 = step( A 3 , b , c, O , 1 , t 1 : echo o f f

title('Stapresp0nsie (kritisch gedempt)') pause, e c h o on plot(t,s3) 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 , XQ31; echo o f f plot( t , r 3

title('1 M o d e (kritisch gedempt)') pause, echo on 1 samengevat echo o f f plot ( title pause plot ( title pause plot ( t, C il, i 2 , i3 1

1

'Impulsresponsies') t , E S I . s 2 , s 3 1 1 'Stapresponsies') t , E rll, r12, r 3 I i

title('2 Bovenkritisch- en 1 kritisch gedempte mode')

pause, e c h o on w h o

Resultaten: b = 0

1

c = 1 0 X Geval I ( b o v e n k r i t i s c h g e d e m p t ) A I = 0 - 4

1

- 8 V l = 1 . 0 0 0 0 -0.5355 DI -0.5359 O -0.1340 1 * 0000 o - 7 . 4 6 4 1 ,>: .?: J '.i '>: -! :: : : i : _ .

A 1 . 1 3

x o l i = 2.0000

- 0 . 5 3 5 9

x o 1 2 = 1 . 0 0 0 0 - 7 . 4 6 4 1

X G e v a l 2 f o n d e r k r i t i s c h g e d e m p t ) A 2 = O - 4 1 - 2 V 2 = -0.2500

-

0.433Oi - 0 . 2 5 0 0 + 0.4330i 1 . 0 0 0 0 1 .O000 O2 =

- 1 .

O000 + 1.7321i 0 o -1.0000

-

1.7321i 1

A 1 . 1 5 1 Gewal 3 ( k r i t i s c h gedempt) A 3 = 0 1 - 4 - 4 V 3 = - 0 . 5 0 0 0 - 0 . 5 0 0 0

1 .oooo

1 .o000 D3 = - 2 a O - 2

x 0 3 = 1 . 0 0 0 0 - 2 . 0 0 0 0

A 1 . 1 7

four variables are: e p s , pi, inf, nan,

...

A l D3 A 2 V I A 3 v 2

o1

v 3 D2 b C i1 i 2 i3 rt 1 r12 i-3 S I s 2 s 3 t U x O 1 1 XOI 2 x 0 3

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.

A

Beschouw daartoe de berekening van e 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)

.

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

PC-flatlab Function-files mreeks.m en eigdia9.m :

t function CE,inorml = mreekst.4)

Z E = Matrix exponential o f 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 = r o u n d ( l o g ( n o r m ( A , ? ) 1 / l o g t 2 ) + 1 . 5 ) ;

i f 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; w h i l e 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)

A 2.4

Resultaten ntet mreeks.m :

A = I O O O O > [ E ? , i n o r m l = m r e e k s f A ) EI = 2.7183 O O O o 2.7183 2.7183 O O O O O 2.7183 U O

inorm = Columns 1 through 7

1 . o 0 0 0

1.0000 O. 2023

Columns 8 through 13

o.

O000

o

* 0000 0 . 0 0 0 0

> log10íinorm)

ans = Columns

1

through 7

O O O 2.7183 O 0.0189 o . O000 O O O O 2.7183 0.0011 o . O000 o . O000 o . OOOO o. O000 O O -0.6941 -1.7244 -2.9685 -4.3656 -5.8792 Columns 8 through 13 -7.4861 -9.1704 -10.9205 -12.7281 -14.5864 -16.4902

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

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

Y ( s ) Y ( s )

Bereken de overdrachtsfunkties

-

en

-

.

*

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.

PC-Matlab Command-file mvsyst.m :

echo o f f , clear, cla

i! 'nvsyst.m' stelt d e matrices M , E , en K o p voor een massa-veersysteem

2 stelt d e matrices A , E , C en R op voor een toestandsbeschrijving

1 berekent d e eigenwaarden van de systeemmatrix

2 berekent d e overdrachtsfunkties

i! maakt Bode- Eon Nyquist-plots van d e overdrachtsfunkties

i!

1 Re invoer wordt opgevraagd. m l = input('rn.i = ' I ; m2 = inputí'rn2 ' I : b l = 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 = l m 1 , 0 ; ~ , m 2 1 K = [kl+k3,-k3;-k3,k2+k33 B = lbl+b3,-b3;-k3,b2ib31 pause A = Izeros(2),eye B = Czeros(2i;inv C = Ceyeí21,zeros D = CzesQsf2)f; pause

A 3.5 CNum1,denll = s s Z t f í A , B , C , D , l ) tNurnZ,den21 = s s P t f ( A , B , C , D , Z ) pause eigenwaarden = eigíA) polen = rootsfden11 pause opnieuw = 1 ; while opnieuw > 0 ,

10 = input('6eef d e logaritme van d e laagste hoekfrekwentie ' 1 : u p = i n p u t ( ' F e e f d e logaritme van d e hoogste hoekfrekwentie ' 1 ;

N = inputí'Aantal plotpunten * I : w = l o g s p a c e í l o , u p , N ) ; i o g 9 0 g ( w , M a g I ) 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

l o g l o g ( w , M a g 2 )

title('Magnitude van overdracht van F Z ' ) pause

semilogxíw,fhase2)

title('face van overdracht van F 2 ' ) pause, grid, pause

axis('square') polaríPhase1,Magl) titleí'ûverdracht van F l ' ) p a u s e , g r i d , p a u s e poi.ar(Phase2,MagZl title('0vesdracht van F 2 ' ) pause, g r i d , pause axis('normal'1

o p n i e u w = input('frekwentiedomein aanpassen ( 1 1 , anders 1 0 ) ' 1 :

end

c l a

e c h o o n w h o

A 3 . 7 Resultaten: ml = 3 m2 = 1 bl = 0 . 3 b 2 = 0 . 8 b3 0.2 k l = 9 k 2 = 5 k 3 = 2 M = 3 O K = 1 1 - 2 B = 0.5000 -0.2000 A = O I -2 7 - 0 n 2000 1 I O000 O O 1 s O000 O O O O 1.0000 - 3 . 6 6 5 7 O . 6667 -0 e 1667 O . 0667 2.0000 -7.0800 o . 2000 - 1 * 0000 8 = 0 O O . 3 3 3 3 O c = 1 O O 1 O O O 1 . o000 O O i i O O

Numí = O

o.

O 0 0 0 O. 3333 O . 3333 2.3333 O

o.

o 0 0 0

o.

O000 O . 0 6 6 7 O. 6667 den1 = 1.0000 1.1667 10.8200 4.5667 24.3333 Num2 = O o . O000 -0.0000 O . 0667 O . 6 6 6 7 O o . O000 I . O000 O . 1667 3.6667 d e n 2 = 2 . 0 0 0 0 1 . 1 6 6 7 1 O . 8200 4 . 5 5 6 7 24.3333 e i g e n w a a r d e n = - 0 . 4 9 6 9 + 2 . 6 6 0 9 i -0.0864 + 1.82031 - 0 . 0 8 6 4 - 1.82031 - 0 . 4 9 6 9 - 2 . 6 6 0 9 1 p o l e n = - 0 . 4 9 6 3 + 2.6603i - 0 . 4 9 6 9

-

2 . 6 6 0 9 i -0.0864

-

1.82031 -0.0864 + 1.8203i

A 3.9 . . . . . < . . . . . . . . . , . , , . . . . . . . , , . . I . . . . . . . . . . . . . . . . . , 1 . . . .

. . . .

A 3 . 1 1

Your v a r i a b l e s a r e : eps, p i , i n f , nan,

...

A B c D fc M Mag1 Mag2 N Numl Num2 P h a s e 1 Phase2 b l b2 b3 m 2 den1 den2 opnieuw elgenwaarden polen k l _UP k 2 W k3 lo m l

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

C = [ 1 O ] 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).

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

Er geldt nu:

figuur 4.2

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.

*

Lz = place( A, B, [ -0.004, -0.006

1

*

xz

= lsim( A-B*Lz, [ O

1

I ' , eye(21, zeros(2,1), wf t ) 2 UZ = -B"Ez*xz'

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

kwadratisch kriterium J = E[ g (t)eQeg(t)

+

u(t)*R*u(t)

1

minimaal is.

Hierin is Q semi positief definiet en R>O.

T

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

feigenwaarden van (A-B*Lo)).

Bereken &,(t) als responsie op w(t) bij terugkoppeling van de toestand

met Lo.

Bereken ook het regelsignaal uo(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.

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

Opdracht 6

Bereken Ks z.d.d. de eigenwaarden van (A-Ks*C) op

-0.500 en -1,000 sek" 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

A 4 . 7

Opdracht 8

Bepaal de optimale terugkoppeling KO van E ( t ) t.a.v. het minimaliseren

van EC fx(t)-s(t) (xtt)-xw(t) 1 I.

Maak hierbij gebruik van Vww en Vw.

Bereken de waarnemerpolen (eigenwaarden van (A-Ko*C)).

Bereken &(t) als responsie op de input van de waarnemer.

Kies als input y,(t) en uo(t).

Plot samen de verplaatsing en de schatting daarvan.

Doe dit ook met de snelheid.

*

lqe

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

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.

PC-Matlab Commandfile:

echo on , clear , cla

1 BATHYSCA.M

2 Optimale regeling van d e duikdiepte van d e bathyscaaf ( R e g e l e n 1 1 ) . 1

2 Berekening van een optimale statische toestands-terugkoppeling en een

1 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 , B = [O;bl, C = [ 1 , 0 3 pause ela

Controllability = ctrbtA,B), Controllabalityrang = rank(Controllabi1ity)

pause cla

Observability = o b s v ( A , C ) , Observabilityrang = rank(0bservability)

pause cla 1 Opdracht 2 t=[0:3003'*3; rand('norma1') w = 0.05*rand(t) ; 1 Systeemruis m / s e c 2

Vww = [ O , O ; O , cov(w)

1

2 Autocorrelatie van d e systeemruis

A 4 . 1 1

echo o f f plot ( t , w )

title('Systeemruis w i t ) Im/secZJ')

pause, echo o n , cla

Z Opdracht 3

eigenwaarden = eigtA)

y = lsim(A, E O , 1 I ' , C , O , w , ti; pause

echo o f f

titlei'Ruisrecponsie zonder regeling') pause, echo on, cla

plotlt,y)

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

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;

u 2 = -Lz*xz' ; us = -Ls*xc' ;

pause

echo o f f

plot( t, xz(:,l) , t , xs(:,ll)

title('Verp1aatsing b i j zwakke en sterke terugkoppeling') pause

plot( t , xzl:,Z) , t , xs(:,21)

title('Sne1heid b i j z w a k k e en sterke terugkoppeling') pause

plot( t ' UZ, t , u s 1

title('Regelsignaa1 b i j zwakke en sterke teTUgkoppeì.ing') pause, echo on, cla

1 O p d r a c h t 5

a = [

1 ,

O ; O , O ] R = 1

pause c l a

2 l m t.t 1 k g

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 = -Lo*xo' ;

pause e c h o o f f

p l o t f t , x o f : , l I 1

title4'Verplaatsing b i j o p t i m a l e t e r u g k o p p e l i n g van d e toestand') p a u s e

plot( t , x o ( : , 2 )

1

title('Sne1heid b i j optimale t e r u g k o p p e l i n g van d e toestand') p a u s s

title('Regelsignaa1 b i j optimale t e r u g k o p p e l i n g van d e toestand') p a u s e , e c h o o n , cla plot( t , u0 1 1 S t a n d a a r d d e v i a t i e s van x en u b i j Z w a k k e , S t e r k e en O p t i m a l e t e r u g k o p p e l i n g SxzSxsSxo = stdl C xz , X S , XO

1 1

S u z S u s S u o = std( [ U Z ' , U S ' , UO' 3

1

pause cla

A 4.13

Z TOESTANDSSCHATTING m.b.v. een WAARNEMER

1. Opdracht 6

K s = ( place( A',C', C -0.500 , -1.000 3 1 1 ' 1 Sterke terugkoppeling B w = E Ks*C , B I ;

xws = lsim( A-Ks*C, Bw, eye(21, meros(2,3ì, X O I Uo'3, t);

pause echo o f f

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 o n , cla

plot( t. xo(:,I) ' t , xws(:,1))

2 Opdracht 7

v = S*rand(t);

z

Meetruis; standaarddeviatie 5 m o p 1 km diepte

v v v = cov(vl Z Autocorrelatie van d e meetruis

Bw = I Ks*C , EL , K S I ;

xwsv = lsimi A-Ks*C, Bw, eye(21, zeros42,4), X o , U O ' , v 1 , t);

pause echo o f $

plot( t , xo(:,l) , t , xwsv(:,.1))

title('Verp1aatsing e n schatting M E T M E E T R U I S ' ) xlabel('SYSTEEM en WAARNEMER NIET gekoppeld') pause

plot( t , xo(:,2) 1 t , xwsw(:,2)1

title('Sne1heid en schatting MET MEETRUIS') xlabel('SYSTEEM en WAARNEMER MIET gekoppeld') pause, echo o n , cla

l Opdracht 6

K O = lqe( A , eye(21, C , Vww, V v v 1 l Optimaal filter

waarnemerpolen = eig( A

-

Ko*C 1

BW = C Ko*C , B , KO

1 ;

xwo = lsim( A-Ko*C, Bw, eye(21, zeros(2,4ì , C xo , uo', v I , t);

pause echo o f f

plot( t , x 0 ( : , 1 ) , t , xwo(:,I))

titlel'Verplaatsing en OPTIMALE schatting MET MEETRUIS') xlabel('SYSTEEí4 en OPTIMALE WAARNEMER NIET gekoppeld') pause

title('Sne1heid en OPTIMALE schatting MET MEETRUIS') xlabel('SYSTEEM en OPTIMALE WAARNEMER NIET gekoppeld') pause, echo o n , c l a

plot( t , X O ( : , 2 ) , t , X W Q ( : , 2 1 ]

2 KOPPELEN VAN DE OPTIMALE WAARNEMER AAN HET TE REGELEN SYSTEEM

l Opdracht 9

Z Matrices van het samengestelde systeem

AZ = [ A , -B*LO ; KO*C, A-B*LO-KO*C .I

82 = [ C O , 1

1'

, zerosi2,I) : zerosi2,l). KO I pause c l a regelpolen,waarnemerpolen gekoppeldepolen = eigfA2I pause c l a

A 4 . 1 5

Z Opdracht 10

x 2 = lsim( A2, 62, eye(41, zeros(4,2), C w , v 1, t);

e c h o o f f

plot( t , x2(:,1) , t , x2(:,3))

title('Verp1aatsing e n schatting')

xlabeli'SYSTEEM en OPTIMALE WAARNEMER gekoppeld') pause

plot( t , x2(:,2) , t , x2(:,4))

title('Sne1heid en schatting')

xlabelf'SYSTEEM en OPTIMALE WAARNEMER gekoppeld') pause p l o t I t i t l e pause plot ( t i t l e pause t , XQ(:,!l , t , X 2 [ : , f ) )

'Verplaatsing b i j terugkoppeling van TOESTAND resp. SCHATTING')

t , xo(:,2) ' t , x2(:,2)1

'Snelheid b i j terugkoppeling van TOESTAND resp. SCHATTING') e c h o on, cla

Z EINDE PROGRAMMA 'bathyscaaf'

Resultaten: a = 0.0100 b = 8.4175E-004 1 O p d r a c h t 1 A = O 0.0001 t . uuoo O B = 1.OE-003

*

O 0 . 8 4 1 8 c =

1

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 . 8 4 1 8 O

A 4 . 1 7 I 1 Opdracht 3 eigenwaarden = 0.0100 -0.0100 7. Opdracht 4 L z = 0 . 1 4 6 5 1 1 . 8 6 0 0 L s = 1 3 . 2 1 6 5 2 6 9 . 4 8 0 2

. . . . . .

A 4.19 1 Opdracht 5 Q = I O O O R = l L O = 1 . 1 2 7 2 5 1 . 1 5 0 7 regelpolen = - 0 . 0 2 1 8 + 0 . 0 1 9 3 1 - 0 . 0 2 1 8

-

0,0193i

SxzSxsSxo = 1 4 4 . 7 5 6 7 O . 6 0 4 8 SUZSUSSUO = 2 4 . 0 5 3 1 1 . 1 7 7 8 3 6 . 6 2 6 2 O . 1 3 1 9 10 - 0 5 2 5 1 7 . 3 5 3 0 O . 2 5 5 2

A 4 . 2 1

1 O p d r a c h t 6 K s = 1 . 5 0 0 0 0 . 5 0 0 1

S Opdracht 7

A 4.23 1 O p d r a c h t 8 K O = 0.1460 0 . 0 1 0 7 w a a r n e m e r p o l e n = -0.0730 + 0.0723i -0.0730

-

0.0723i

2 O p d r a c h t 9 0 0 0.0001 O - 0 . 0 0 0 3 -0.0436 0.1460 O -0.1460 1.0000 0.0107 O -0.0115 -0.0436 A2 = O 1.0000 O O O O 0.1460 O 0.0107 82 = 1 .o000 r e g e l p o l e n = -0.0218 + 0.0193i -0.0218 - 0.0133i waarnemerpolen = -0.0730 + 0.0723i -0.0730

-

0.07231 g e k o p p e l d e p o l e n = e i 9 (A2 1 g e k o p p e l d e p o l e n = -0.0730 + 0.0723i -0.0730

-

0.0723i -0.0218 + 0.0193i -0.8218

-

0.0193i X Opdrach t 10 .-. i. i i.

-

.,.: j

Y o u r variables are: e p s , pi, inf, nan,

...

A A 2 8 8 2 8 W C Controllability Controllabilityrang KO K S Lo L s L Z Observability Observabilityrang Q R suzsussuo sxZsxssx0 v v v VWW a b eigenwaarden gekoppeldepolen regelpolen t U 0 us UZ V W waarnemerpolen x2 XO x s xwo xws xws v XZ Y