Het toepassen van semidefinite programmering bij het oplossen van regeltechnische problemen

Het toepassen van semidefinite programmering bij het

oplossen van regeltechnische problemen

Citation for published version (APA):

van de Laar, J. H. A. (1996). Het toepassen van semidefinite programmering bij het oplossen van regeltechnische problemen. (DCT rapporten; Vol. 1996.113). Technische Universiteit Eindhoven.

Document status and date: Gepubliceerd: 01/01/1996

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

TECHNISCHE UNIVERSITEIT EINDHOVEN

Faculteit Werktuigbouwkunde

Vakgroep Fundamentele Werktuigkunde (WFW)

Sectie Systeem-Regeltechniek

Stageverslag

:

Begeleider : &ir.

A.G. de Jager

Het toepassen van Semidepnite

Programmering bìj het oplossepz

van regeltechnische problemen

J.H.A.

van de Laar

augustus 1996

Inhoud

Samenvatting 111

1 Hoofdstuk 1 : Inleiding

Hoofdstuk 2 : Overzicht van de geschiedenis van de LMI’s en SP 2

tj 2.1 : Geschiedenis van LMl’s 2

tj 2.2 : Geschiedenis van de SP methode 3 5 Hoofdstuk 3 : Semidefinite Programmering

5

3

3.1 : Introductie

tj 3.2 : Dualiteit

tj 3.2.1 : Inleiding

tj 3.2.2 : Het duale SP

tj 3.2.3 : Omschrijving van het primair-duale probleem

tj 3.3 : Het centrale pad

tj 3.3.1 : Inleiding

tj 3.3.2 : Grensfunctie voor een LMI

3

3.3.3 : Analytisch centrum van een LMI

tj 3.3.4 : Berekening van het analytisch centrum

tj 3.3.5 : Parametrisatie van het centrale pad

9 9 9 10 10 11 12 12 12 14 14

tj 3.4 : Methoden voor het oplossen van SP’s

3

3.4.1 : Inleiding

tj 3.4.2 : Primair-duaal potentiaal reductie methoden

tj 3.4.3 : Lange-stap pad-volgende methoden

tj 3.4.4 : Plane search

tj 3.5 : Phase I en Phase I-Phase I1 methoden 15 17 Hoofdstuk 4 : Software tj 4.1 : Inleiding 17 tj 4.2 : Het sp-pakket tj 4.2.1 : sp.m tj 4.2.2 : Convergentie criteria tj 4.2.3 : Opmerkingen 17 17 18 19 tj 4.3 : LMItool tj 4.3.1 : De functie 1mitool.m

tj 4.3.2 : De werking van LMItool

19 20 20

Hoofdstuk 5 : Het werken met de LMI methode

0

5.1 : Inleiding

5

5.2 : Het LQR probleem

tj 5.3 : Het inverse probleem van de optimale regelaar Hoofdstuk 6 : Conclusies en aanbevelingen

Literatuurlijst Afkortingen Notaties Bij lagen 22 22 26 31 33

Samenvatting

Een nieuwe methode voor het oplossen van regeltechnische problemen is het toepassen van Semidefinite Programmering. Bij deze SP methode wordt het probleem geschreven als een SP. Hieronder staat zo’n SP weergegeven

minimaliseer C T X

onder de voorwaarde F ( x ) 2 O

IE een SP wor& eefi ~pei-ratiedoe~

(~5)

gedefinieerd en er worden eïkele eisen gesteld

waar de variabelen aan moeten voldoen. Met deze eisen wordt een gebied afgebakend waarbinnen de waarden voor de variabelen moeten liggen. Als alle variabelen binnen het afgebakende gebied liggen is het probleem feasible. Bij (1) hoort het hieronder gegeven duale probleem

maximaliseer - Tr F,Z

onder de voorwaarde Tr E;]Z = ei, i = 1,. . ., m 2 2 0

De primair-duale potentiaal reductie methode beschrijft een algoritme waarmee SPY s kunnen worden opgelost. Bij deze methode wordt een potentiaal functie p gereduceerd met een vast bedrag per stap. Bij deze methode blijfi de iteratie feasible en convergeert de iteratie naar een optimum.

Een ander algoritme voor het oplossen van SP’s wordt beschreven door de lange-stap pad-volgende methode. Deze methode bestaat uit twee verschillende stappen. In de eerste stap, de centralisatie stap, wordt de afwijking van centraliteit tot nul geredu- ceerd, terwijl het dualiteitsverschil constant wordt gehouden. In de tweede stap, de dualiteit sverschil-reductie stap, wordt het dualiteitsverschil geminimaliseerd in een vlak waarvoor geldt dat de afivijking van centraliteit kleiner is dan een van tevoren vastge- stelde waarde.

De beginwaarden dienen strikt primair en duaal feasible te zijn. In de meeste gevallen wordt aan deze eis voldaan. Zijn de beginwaarden echter niet strikt primair en duaal feasible dan kan het SP worden aangepast, zodat er toch een oplossing gevonden kan

worden. Een van de methoden waar dit mee gedaan kan worden is de big-Mmethode. Bij de big-Mmethode wordt een extra eis toegevoegd aan het SP. Het spoor van F(x)

(1) mag niet groter zijn dan de waarde M. Als M groot genoeg gekozen wordt dan heeft deze extra voorwaarde geen invloed en wordt in wezen het oorspronkelijk probleem (1) opgelost.

De software die gebruikt kan worden bij het oplossen van SP’s bestaat uit het sp- pakket en LMItool, een interface voor het sp-pakket. Het sp-pakket bestaat uit

verschillende functies waarmee SP’s zoals (1) kunnen worden opgelost. Het program-

ma geeft een oplossing als een van de (vijf) convergentie criteria bereikt wordt. Voor deze convergentie criteria zijn het aantal iteraties, het dualiteitsverschil en de target value van belang.

LMItool fungeert als een intefiace voor het sp-pakket. Voordat een van de functies uit het sp-pakket kan worden aangeroepen moeten er veel variabelen worden gedefi- nieerd. LMItool neemt dit werk voor een groot gedeelte over. Met LMItool worden twee functies gegenereert. In een van deze functies, de evaluatie functie, moet het op te lossen SP gedefinieerd worden aan de hand van de in te voeren LMI’s en het operatiedoel.

Bij het oplossen van het LQR probleem met de SP methode volgen dezelfde waarden

als bij het oplossen van het LQR probleem met een klassieke methode. Het oplossen

vm het LQR probleem met de SP mthode is echter &et voordeliger ten opzichte van

het oplossen met de klassieke methode. Bij systemen van dimensie 4 is het verschil in benodigde oplostijd al behoorlijk. Voor grotere systemen zal dit verschil alleen maar toenemen.

Een voordeel dat de SP methode wel heeft is dat regeltechnische problemen kunnen

worden opgelost, welke met klassieke methoden niet oplosbaar zijn. Voor het oplossen van het inverse probleem van de optimale regelaar is geen klassieke methode bruik- baar. Na controle van de gevonden waarden blijkt met de SP methode een goede

oplossing gevonden te worden.

Op dit gebied is de SP methode een aanvulling van het assortiment oplosmethoden voor regeltechnische problemen.

Als een ïegeltechisch probleein opgelost kan woïdeil met zowel eei2 klassieke

methode als met de SP methode, dan heeft de klassieke methode de voorkeur boven de

1 : Inleiding

Dat echt alle wegen naar Rome leiden is erg onwaarschijnlijk, maar een heleboel toch zeker. Voor de regeltechniek geldt in ieder geval dat er vele wegen naar Rome leiden. Een van de nieuwe wegen is het oplossen van regelproblemen door gebruik te maken van de SP methode.

Doel van dit onderzoek was kijken naar de mogelijkheden die de SP methode heeR in het oplossen van regelproblemen. Hierbij is gekeken naar toepassingen in zowel problemen &e met H29sieke ~ ~ e t h ~ c l e f i oglosbaar zijn, als p r ~ h ! e ~ e ~ w a r nog geen

algoritme voor bekend is.

Met de SP methode is het mogelijk om vraagstukken met de computer op te lossen waar voorheen nog geen efficiënt algoritme voor bestond. Daarnaast kunnen ook vraagstukken opgelost worden waar al wel een algoritme voor bestaat.

Bij de SP methode wordt het probleem omschreven met LMI’s.

Na een korte beschrijving van de geschiedenis van de LM’s en van de SP methode in hoofdstuk 2, zal in hoofdstuk 3 de technische achtergrond van de SP methode worden

belicht. Hierin komt aan de orde welk algoritme gebruikt wordt en zal worden aange- geven hoe dit algoritme wedd. in hoofdstuk 4 komt de sofiwzrs aan bod waarmee

SP’s opgelost kunnen worden. Hierin worden de belangrijkste parameters van het

hoofdprogramma sp.m aangestipt. Verder wordt het programma LMItool belicht, LMItool is een gebruiksvriendelijke interface voor het programma sp.m.

Na deze theoretische achtergrond van het verhaal worden in hoofdstuk 5 twee voor- beelden als SP uitgewerkt. In § 5.1 wordt het LQR probleem bekeken om de SP methode te vergelijken met een klassieke methode. Vervolgens wordt in

5

5.2 het inverse probleem van de optimale regelaar belicht. Voor het oplossen van dit probleem met de computer zijn nog geen geschikte, klassieke algoritmen bekend.

Naar aanleiding van de resultaten uit deze voorbeelden volgen tenslotte in hoofdstuk 6 conclusies en aanbevelingen.

2 : Ovenicht van de geschiedenis van de LMI’s en de SP methode 2.1 : Geschiedenis van LMI’s

De geschiedenis van de LMI’s begint rond 1890, bij de introductie van de Lyapunov theorie. Lyapunov liet zien dat de differentiaalvergelijking

d dt

- x ( t ) = A x ( t )

alleen stabiel is als voldaan wordt aan de eis dat er een positief-definite matrix P

bestaat zodat

De eis P > O, ATP

+

PA < O wordt een Lyapunov ongelijkheid in

P

genoemd en is een speciale vorm van een LMI.

Lyapunov liet zien dat deze LMI expliciet kon worden opgelost door een Q = Q‘ > O

te nemen en vervolgens de vergelijking ATP

+

PA = -Q voor de matrix P op te lossen,

P

is dan gegarandeerd positief-definiet als het systeem (2.1) stabiel is. De eerste LMI die gebruikt werd om de stabiliteit van een dynamisch systeem te analyseren was de Lyapunov ongelijkheid (2.2), welke analytisch kan worden opgelost door een verzameling lineaire vergelijkingen op te lossen.

In de jaren ‘40 paste Lur’e, Postnikov en anderen in de Sowjet Unie Lyapunov’s methoden toe op specifieke, praktische problemen in de regeltechniek. Ook al werden door hen niet expliciet LMI’s gevormd, toch hadden hun stabiliteitscriteria wel de vorm van LMT’s. Lur’e en anderen waren de eersten die Lyapunov’s methoden

toepasten op praktische regelproblemen. De LMI’s die hier uit volgden werden, met de

hand, analytisch opgelost.

Een volgende belangrijke stap werd in de jaren ‘60 gezet. Yakubovich, Popov, Kalman en andere onderzoekers slaagden erin de oplossing van de LMI’s, welke voortkwamen uit het probleem van Lur’e, te reduceren tot simpele grafische criteria door middel van het positive-real

(PR)

lemma [Bi]. Dit resulteerde in het Popov-criterium, Tsypkin

criterium en vele variaties.

Men liet zien hoe een bepaalde groep van LMI’s kan worden opgelost door middel van grafische methoden.

Het positive-real

(PR)

lemma en uitbreidingen werden intensief bestudeerd in de tweede helR van de jaren ‘60, Men ontdekte dat deze in verband stonden met de ideeën over passiviteit, de small-gain criteria en kwadratisch optimale regelingen.

Rond 1970 was bekent dat de LMI, welke voorkomt in het PR lemma, niet alleen via grafische weg opgelost kon worden, maar ook door het oplossen van een algebraïsche ficcati vergelijhg

(Ml.

A ~ P + P A + Q

m+cT

B T P + C R

kan worden opgelost door te kijken naar de symmetrische oplossingen van de ARE

Deze oplossingen kunnen worden gevonden via een eigenwaardendecompositie van de bijbehorecde Hami!ton-mtrix.

Rond 1970 waren dus diverse methoden bekent voor het oplossen van bepaalde speciale typen LMI’s : direkte (bij kleine systemen), grafische methoden en door het oplossen van Lyapunov of Riccati vergelijkingen.

De volgende grote vooruitgang was de constatering dat de LMI’s die voorkomen in systeem- en regeltechniek kunnen worden geformuleerd als convex O ~ ~ i m ~ Z i s e ~ i n g s - problemen (voor de definitie van de term convex zie

[CAM],

[ROC]) die geschikt zijn

voor het oplossen met computer. Met deze constatering is het nu ook mogelijk LNII’s op te lossen die geen analytische oplossing kennen of waarbij een analytische oplossing moeilijk te vinden is.

Pyainitsky eíì Sksrodinsky [PSI warez in de jarez ‘80 de eerste wetefisctlappers die dit helder en duidelijk bewezen. Zij waren de eersten die een Lyapunov functie als convex optimaliseringsprobleem formuleerden, en daarna een algoritme toepaste waarmee altijd een oplossing voor het optimaliseringsprobleem verkregen werd.

Van groot praktisch belang is de ontwikkeling van een sterke en effectieve interior- point methode voor het oplossen van LMI’s die voorkomen in de systeem- en regeltechniek. In 1984 introduceerde Karmarkar een nieuw lineair programmeer algoritme die lineaire programma’s in polynominale tijd oplost en die in de praktijk zeer efficiënt bleek te zijn. In de tijd hierna was het onderzoek toegespitst op algoritmen voor lineaire en (convex) kwadratische programma’s.

In 1988 ontwikkelden Nesterov en Nemirovsky interior-point methoden die toepasbaar zijn bij het oplossen van convex optimaliseringsproblemen.

Ook al is er nog veel te doen op dit gebied, toch zijn al enkele interior-point algoritmen voor LMI problemen geïmplementeerd en getest op bepaalde LMI’s die voorkomen in de regeltechniek. Hierbij zijn ze zeer efficiënt gebleken.

§ 2.2 : Geschiedenis van de SP methode

[BF] was in 1963 het eerste artikel dat de theoretische eigenschappen van de

semidefinite programmering behandelde. Veel wetenschappers hebben in de jaren ‘70

en ‘80 gewerkt aan het minimaliseren van de maximale eigenwaarde van een

symmetrische matrix welke geschreven kan worden als een SP zie [GT,OWl,OWZ]. De interior-point methode is een nog relatief jonge methode. Interior-point methoden voor LP’s werden in 1984 geïntroduceerd door Karmarkar. Dit algoritme, en de methoden die later

ontwikkeld

werden, combineren een lage, polynominale worst-case

complexiteit met uitstekend gedrag in de praktijk. Interior-point methoden werden daarna uitgebreid om convex kwadratisch programmering aan te kunnen.

Een belangrijke doorbraak werd bewerkstelligd door Nesterov en Nemirovslq in 1988 [NNI. Hierin lieten zij zien dat interior-point methoden voor lineair programming, in principe, gegeneraliseerd konden worden naar alle convexe optimali-seringsproblemen. Ook Kamath en Karmarkar [KK] lieten deze mogelijkheid zien. SP's zijn een

belangrijke groep van convexe optimaliserings-problemen waarvoor interior-point methoden bruikbaar zijn.

3 : Semidefinite Programmering

In dit hoofdstuk wordt de technische achtergrond van de SP methode belicht. Er komen verschillende termen aan bod die voorkomen in het algoritme dat wordt gebruikt bij het oplossen van SP's. Dit hoofdstuk is geschreven aan de hand van [BZ].

Voor meer informatie over de technische achtergrond van de SP methode wordt dan ook verwezen naar [B2].

5

3.1 : Introductie

We bekijken het probleem van het minimaliseren van een lineaire functie van een variabele

x

E R" die voldoet aan de LMI F(x) 2 O :

T minimali seer c x onder de voorwaarde F( x) 2 O met M F ( x ) F,

+

xJ;; i=l xi E X , i = l ) "", m

en cTx het primaire operatiedoel.

De vector c E R" en m + l symmetrische matrices Fa,

...,

F, E R"" vormen de data van het probleem. F(x) 2 O is een M en 2 wil zeggen dat F(x) positief-semidefiniet is. Daarom wordt (3.1) een SP genoemd.

Omdat het doel en de voorwaarde convex (zie [CAM], [ROC]) zijn wordt een SP zoals (3.1) een convex optimaliseringsprobleem genoemd.

Het voorbeeld in figuur 1 geeft een idee van de kenmerken van een SP met

x

E

R2

en

F, E

R7x7.

t

i

F(x)r

O

I

De doorgetrokken lijn is de grens van het feasible gebied. Het feasible gebied wordt gevormd door de grenscurve en het gebied dat er door wordt ingesloten. Het SP probleem komt er eigenlijk op neer om zo ver mogelijk in de richting -c te lopen zonder het feasible gebied te verlaten. Bij dit SP is er één optimaal punt, xopf. In dit voorbeeld komen enkele algemene kenmerken van SP’s naar voren. Zoals al eerder vermeld is het SP convex [ C M R O C ] . Het optimale punt ligt op de grens van

het feasible gebied, F(xopJ is singulier, in het algemeen ligt (bij een feasible probleem)

het optimaal punt altijd op de grens van het feasible gebied.

In dit voorbeeld is de grens van het feasible gebied niet glad . Bij verreweg de meeste gevallen bestaat de geris

u:iL

stuksgewijs 21gebrgsche oppervlakken.

De SP methode kan worden beschouwd als een uitbreiding van lineair programmeren. Het LP

minimaliseer cTx

onder de voorwaarde

Ax

+

b 2 O (3.3)

kan worden geschreven als een SP met F(x) = diag@x+b).

Het is dan ook niet verwonderlijk dat de theorie van de SP methode overeenkomsten

vertoont met de theorie van lineaire programmering.

Er zijn echter belangrijke verschillen. Zo zijn de resultaten met betrekking tot de dualiteit zijn voor SP’s zwakker dan voor LP’S.

Het volgende voorbeeld van een niet-lineair optimaliseringsprobleem kan worden geschreven als een SP, maar niet als een LP.

minimaliseer ( C T x y

d T x

onder de voorwaarde Ax

+

b 2 O

aangenomen wordt dat cfx > O als

Ax

+

b 2 O.

Allereerst wordt een hulpvariabele t geintroduceerd die dient als bovengrens van het primaire operatiedoel : mimimaliseer t onder de voorwaarden

Ax

+

b 2 O I t ( C ’ x y dTX (3.5)

In deze vorm is het primaire operatiedoel een lineaire funktie van x en f. Het niet- lineaire operatiedoel van (3 -4) is nu een niet-lineaire voorwaarde in (3.5). Door middel van Schar complements [B 1, p. 281 worden de voorwaarden in (3.5) omgeschreven

minimaliseer t onder de voorwaarde cTx dTx diag(Ax+b) O O O

Het niet-lineaire probleem (3.4) is nu omgeschreven naar een SP 43.6).

In dit voorbeeld zijn twee trucs gebruikt die vaak toegepast worden. Het schrijven van meervoudige LMI’s als een blokdiagonale matrix ongelijkheid en het gebruiken van

Schw compiemenis am eeïì ï i e t - h e ~ r e c m v e ~ e vco~s~2ardei. te schijven als eer? LMI.

Het herkennen van Schur complements in niet-lineaire voorwaarden is vaak de sleutel in het omschrijven van niet-lineaire convexe optimaliseringsproblemen naar SP’s.

Het feit dat positief-semidefinite voorwaarden direkt voorkomen bij vele toepassingen en de mogelijkheid om (vele) convexe optimaliseringsproblemen te schrijven als SP’s

zijn enkele redenen die genoemd worden om de SP methode te bestuderen. De SIP methode biedt dus de mogelijkheid om de eigenschappen van en algoritmen voor een reeks van convexe optimaliseringsproblemen te bestuderen. De belangrijkste reden om de SP methode te bestuderen is dat SP’s zowel in theorie als in praktijk effectief opgelost kunnen worden.

De mate waarin de SP methode een LLconcurrent” is voor klassieke oplosmethoden komt terug in de uitgewerkte voorbeelden.

In [B2] wordt bij het oplossen van SP’s het gebruik van interior-point methoden geprevereerd boven andere bruikbare methoden vanwege de praktische effectiviteit, theoretische effectiviteit en de mogelijkheid om de probleemstructuur te exploiteren.

3.2 : Dualiteit

5

3.2.1 : Inleiding

Een deel van het algoritme waarmee SP’s worden opgelost bestaat uit het

minimaliseren van het dualiteitsverschil. Het dualiteitsverschil is aaiankelijk van het primaire operatiedoel en het duale operatiedoel. Het primaire probleem is reeds in het voorgaande behandelt. In

0

3.2.2 zal het duale probleem worden belicht. Minimalisatie

van het dualiteitsverschil komt aan de orde in

0

3.2.3, waar wordt gekeken naar het oplossen van het primair-duale optimalisatie probleem.

§ 3.2.2 : Het duale SP

Het duale probleem van het SP (3.1) is

maximali seer

-

Tr F3Z

onder de voorwaarde Tr F,Z = ci i = 1,. . . m 2 2 0

met 2 als variabele.

Dit duale probleem is ook een SP, het kan in dezelfde vorm worden gegoten als het primaire probleem. In het algemeen kan het duale probleem vereenvoudigt worden als de matrices F, een bepaalde structuur hebben. Als de matrix F(x) een blokdiagonaal is mag worden aangenomen dat ook de duale variabele Z de structuur van een blokdiago- naal heeft. De speciale structuur leidt tot een vereenvoudiging van het duale probleem. De belangrijkste eigenschap van het duale probleem is dat het de optimale waarde van het primaire SP bepaald. De waarde van het duale operatiedoel van een duaal feasible

punt Z is kleiner of gelijk aan de waarde van het primaire operatiedoel van een primair feasible punt. Het verschil tussen deze twee wordt het dualiteitsverschil q genoemd

q 2 cTx

+

Tr F,Z = Tr F(x)Z _(3.8) Laatp* de optimale waarde van het SP (3.1) zijn p *

d“

de optimale waarde van het SP (3.7) d* 2 sup{-Tr F,Z

I

Z = ZT 2 O, T r q Z = c,} . Stelling 1 [B2, $3. i] :

inf(cTx

I

F ( X ) 2 O }en

Er geldt p * = d als aan een van

de

volgende voorwaarden volähan is

I . Het primaire probleem (3.1) is strikfeasible, er is een x met F(x) 2. Het &ale probleem (3.7) is strik-feasible, er is een Z met Z = ZT

û

O, Tr F,Z= Ci Stelling 1 geeft voorwaarden voor optimalisering van het SP (3. i), met strikt primaire en duale feasibility : x is optimaal alleen als er een Z bestaat zodat geldt

F ( x ) 2 O

Z 2 O, Tr

42

= ei, i = 1,. . . . , m

ZF(x) = O

(3.9)

In [B2, $3.11 zijn enkele voorbeelden met bovenstaande stelling opgenomen.

5

3.2.3 : Omschriiving; van het primair-duale probleem

Voor de meeste SP’s die worden bekeken (welke worden opgelost met interior-point

methoden) geldt één van de voorwaarden van stelling 1 , p*=

d“

(zo niet, dan is het

mogelijk het probleem te reduceren tot een probleem waarvoor een van de voor- waarden wel geldt, zie fB1, $2.51).

Primair-duaal methoden voor SP’s genereren een reeks van primaire en duale feasible

punten

x”)

en .Z@’ (met k = O,l, ... ).

x”)

wordt geïnterpreteerd als een suboptimaal punt welke de bovengrens geeft : p* 4 cTx”) en

Zn”‘

bepaald de ondergrens : p * 2 -Tr F@).

De suboptimaliteit van x” kan worden begrensd in termen van het dualiteitsverschil Als stop-criterium wordt

q”) : c ~ x ( ~ ) - p* I _{q(k) cTxCk) + TrFoZ(k)}

c T d k )

+

Tr F,Z‘’’ I E (3.10) genomen met E > O , een opgegeven tolerantie. Hiermee wordt bij de uitgang een E

-

In I132, $3.21 staat een voorbeeld waarmee dit idee wordt geïllustreerd.

In dit voorbeeld wordt het onderscheid duidelijk tussen convergentie naar een gegeven tolerantie en het garanderen van convergentie tot een gegeven tolerantie.

Primair-duaal methoden kunnen worden gezien als het oplossen van het primair-duale optimaliseringsprobleem minimali seer c T x

+

TrF,Z onder de voorwaarde F ( x ) 2 O, 2 2 O - 7 i r b i Z = ci, i = i, ....? m (3.11)

Hier wordt het dualiteitsverschil cTx

+

Tr F&’ geminimaliseerd over alle primaire en duale feasible punten, de optimale waarde is nul. Het dualiteitsverschil is een lineaire functie van x en Z en is daarom te schrijven als een SP in x en Z. Het voordeel van het bekijken van primair-duale optimaliseringsprobleem (3.1 i) ligt in het feit dat in elke stap duale informatie wordt gebruikt om een nieuwe waarde voor x”l te vinden, en

andersom.

jj 3.3 : Het centrale pad

Vanaf nu wordt uitgegaan van strikt primaire en duale feasibility bovendien worden de matrices F,, i= 1, ..., m lineair onafhankelijk verondersteld.

jj 3.3.1 : Inleiding

Een ander deel van het algoritme waar SP’s mee worden opgelost behelst de

minimalisatie van de afivijking van het centrale pad. Het primaire centrale pad is opgebouwd uit de punten x*(d van het analytisch centrum. Waar yde waarde van het primaire operatiedoel is , y E [ p * ,

p l .

Het analytisch centrum is de waarde x* waarvoor

de grensfunctie geminimaliseerd wordt.

Via

F(x*(d)-’/iz

is het duale centrale pad indirect ook parametriseerbaar met y . Met het oog op het beschrijven van primair-duale algoritmen, komt het beter uit om het centrale pad te parametriseren met het dualiteitsverschil.

jj 3.3.2 : Grensfunctie voor een LMI

De functie

log det F(x)-’ if F ( x ) > O sa0 anders

(3.12)

is een grensfunctie voor X

eindig als

p;O

> O en oneindig ais x de grens vanX nadert.

In figuur 2 zijn de contouren van de grensfunctie van het SP uit figuur 1 te zien.

Figuur 2

De functie is vlak in het midden van het feasible gebied en neemt bij de grens snel toe.

5

3.3.3 : Het analytisch centrum van een LMI

X is begrensd. Daar

Kx)

convex [CAM,ROC] is heeft het een uniek minimum

x* 2 arg min &O) . x* wordt het analytisch centrum van de LMI F(x) 2 O genoemd. Waarbij x* afhankelijk is van de LMI en niet van het oplosgebied X . Wanneer namelijk overtoIlige beperkingen opgelegd worden aan het SP uit figuur 2 zal het analytisch centrum naar een andere positie verschuiven.

x* wordt gekarakteriseerd door TrF(x*)-'

F]

= O, i = 1,. . . , rn

dit komt overeen met de gelijkheidsvoonvaarde Tr E Z = c, , i = I,. . . , m welke voor- komt in het duale SP (3.7). Dat er een verband bestaat tussen duale feasibility en analytische centra zal later duidelijk worden.

3

3.3.4 : Berekening van het analvtisch centrum

De berekening van het analytisch centrum welke hier kort zal worden aangestipt is een voorbode van de primair-duaal interior-point algoritmen die later aan de orde komen. Gegeven een strikt feasible x is het algoritme

__

m

1. bereken de Newton richting

@,

ckN =

3. pasxaan: x : = x + p h N

totdat x, met de gegeven nauwkeurigheid, niet meer afneemt.

De norm

11.11~

is de Frobenius-norm [B2,54.3], stap twee kan worden gedaan met standaard methoden zoals bisectie.

In

m,

52.21 wordt een uitgebreide analyse van de snelheid van convergentie bij dit algoritme gemaakt.

Stelling 2

-Met x”) de wuarde van x in bovenstaand algoritme na k iteraties, en O voor E I . Geldt (3.13) 1 k 2 3.26(@(~“’) -

#(x*))

+

log, log,(-) E Voor de ~ ~ k e u r i g h e i d g e l d t h r t #.x”)j-@(*)j C E

Het is opvallend te zien dat (3.13) onafhankelijk van de grootte van het probleem is. 3.3.5 : Parametrisatie van het centrale pad

Gegeven de LMi

F ( x ) > O

T

c x = y (3.14)

met p * < y <

p

2 sup{cTx

I

F ( x ) > O } . Onder de aanname dat SP (3.1) strikt primair en duaal feasible is wordt het analytisch centrum van (3.14) gegeven door

x*( y ) arg min log det F(x)-’

onder de voorwaarde F ( x ) > O (3.15)

T

c x = y

De lijn beschreven door x*@ wordt het centrde pad genoemd.

Hier is het centrale pad geparametriseerd door de primaire operatiedoel

Bijna alle interior-point methoden benaderen het optimale punt door het centrale pad te volgen. Voor strikt feasible x kan de ahijking van het centrale pad ~ ( x ) worden geschreven als

~ ( x ) log det F(x)-’ - log det F(x*(cTx))-’ (3.16)

~ ( x ) is het verschil tussen de waarde van de grensfunctie bij x en het minimum van de

grenshunctie gezien over alle punten met dezelfde waarde voor de kostfunctie als x.

Bij het beschrijven van primair-duale algoritmen is het parametriseren van het primaire centrale pad en het duale centrale pad met het dualiteitsverschil het meest geschikt

(x*( q), Z*

(v))

f arg min - log det F ( x )

-

log det Z onder de voorwaarde F ( x ) > O, Z > O Trl;TZ= ci, i = l,...,m c T x

+

Tr F,Z = 7 (3.17) voor q

O

Het ver9chi! (3.18)

v ( x , Z ) -log detF(x)Z

+

log det F(x*(v))Z*(q)

= -log det F(x)Z

+

nlog Tr F ( x ) Z - nlog n

is een maat voor de afwijking van

(x,z)

van centraliteit : ~ ( x , 2 ) is een bovengrens voor de inspanning die nodig is om

(x,z)

te centreren.

y(x, 2 ) is alleen afhankelijk van de eigenwaarden &,..,&van F(x)Z n

(3.19)

i=l

3.4 : Methoden voor het oplossen van SP’s 3.4.1 : Inleiding

In het voorgaande zijn enkele begrippen behandeld die voorkomen in de algoritmen van de verschillende methoden die gebruikt kunnen worden voor het oplossen van

SP’s. In deze paragraaf worden deze methoden nader belicht. Een groep van deze

methoden wordt gevormd door de potentiaal reductie methoden. In

3

3.4.2 wordt een algemene omschrijving gegeven van de potentiaal reductie methode. Voor de

verschillende versies van de potentiaal reductie methode wordt verwezen naar

[BZ,

De lange-stap pad-volgende methode kan ook gebruikt worden voor het oplossen van SP’s. Deze methode wordt belicht in

6

3.4.3.

55.2-5.51.

6

3.4.2 : Primair-duaal potentiaal reductie methoden

Hieronder wordt een algemene beschrijving gegeven van de potentiaal reductie methode.

Potentiaal reductie methoden zijn gebaseerd op de potentiaal functie

q ( x , 2) 2 v&log(Tr F ( x ) Z )

+

~ ( x , 2)

(3.20) (n

+

v&)log(Tr F ( x ) Z ) - log detF(x) - log det Z - nlog n

hierin wordt het dualiteitsverschil van het paar x, 2 gecombineerd met de afwijking van centraliteit van het paar x, Z. De constante v > 1 bepaald het gewicht van de term met

het dualiteitsverschil en die met de afijking van centraliteit.

Potentiaal reductie methoden beginnen met een strikt feasible paar

do)

en

Z'O)

waarna p

wordt gereduceerd met minstens een vast bedrag per stap

waarin 6 een positieve constante is.

Diî heeft als gevolg dat de iteratie feasible SlgR en naar een o p t L i ~ m convergeert. Stelling 3 [B2, $5.11 :

Aangenomen wordt dat (3.21) geldt voor S

met O < E I

d m

geldt voor

1

O welke niet aflangt van n eng

v&log(-)

+

y(x'o',z(o))

E

k > dat Tr F(x"/)Z@) rTr F(X(~))Z(~)

s

In het algemeen ziet het potentiaal reductie algoritme er als volgt uit :

Gegeven strikt feasible x en Z.

Herhaal

1. vind een geschikte richting &, en een geschikte duale feasible richting

a.

2. vind eenp, q E R zodat q$x+pSx,Z+qSZ). 3. pas x aan : x: = x+p& en Z : = Z + q a . totdat dualiteitsverschil 7 I t:

De primaire richting & volgt uit

De duale richting

SZ

volgt uit

m -S(D

+

C

aiST<S)ST i=l (3.22) (3.23) (3.24)

D =

DT

en S hangen af van het specifieke probleem en veranderen bij iedere iteratie

Hierboven is het algemene idee weergegeven van potentiaal reductie methoden. In [B2, $5.2-5.51 worden verschillende potentiaal reductie methoden behandelt met verschillende achterliggende gedachîen.

3

3.4.3 : Lange-stap pad-volgende methoden

Bij een potentiaal reductie methode wordt bij iedere stap de potentiaal functie gereduceerd . De potentiaal functie bestaat uit termen met het dualiteitsverschil en afwijking van centraliteit. Het doel is dualiteitsverschil te reduceren terwijl de afwijking van centraliteit gelijk wordt gehouden. De parameter v bepaald inwezen de grootte van de toegestane at'wijjking van centraliteit.

In lange-stap pad-volgende methoden wordt hetzelfde bereikt als bij de potentiaal reductie methoden. Het verschil is alleen dat dit gebeurt in twee verschillende stappen. De eeïste s t q is eeiì cc~âí-a!isatie st2p, hierin wordt de afixijhg vm centrabteit tot nul gereduceerd terwijl het dualiteitsverschil constant wordt gehouden. In de tweede stap, dualiteitsverschil-ïeductie stap, wordt het dualiteitsverschil geminimaliseerd in een vlak waarvoor geldt dat de afivijking van centraliteit

(v)

kleiner is dan de maximale waarde. Het resultaat is een ahame van het dualiteitsverschil waarbij de ahijking van

centraliteit beperkt wordt .

Het algoritme ziet er als volgt uit :

Gegeven strikt feasible x en Z met d x , Z ) 5 ymm

Herhaal

1. Centralisatie stap : centreer x, Z terwijl het dualiteitsverschil constant blijft

2. Dualiteitsverschil reductie stap :

a) Bereken richtingen

Wd,

&Yd loodrecht op het centrale pad

b) Vindp, q E R zodat Tr F ( x + p W ~ ( Z + q P ~ geminimaliseerd wordt

c) pas x,

z

aan : x:= x + p s d ,

z

:=

z + ~ s z P ' ~ .

met ~V(X+PH~,Z+~S~ 5 ymm. totdat dualitietsverschil q 5 E.

In @32, $5.6.31 is een voorbeeld van deze methode te vinden.

9

3.4.4 : Plane search

Zowel bij de potentiaal reductie methoden als bij de lange-stap pad-volgende methode wordt een primaire richting & en een duale richting

SZ

geselecteerd. De lengte van de stappen in de desbetreffende richtingen wordt vastgesteld met een plane search. Het zal blijken dat de inspanning die plane search gedaan moet worden sterk geredu- ceerd kan worden door de matrices eerst te diagonaliseren. Met deze aanpassing is de inspanning nodig voor de plane search een verwaarloosbaar kleine ten opzichte van de totale inspanning.

Voor de hierboven beschreven lange-stap pad-volgende methode is het plane search probleem :

(3.25) minimaliseer pcT

Loc

+

qTrF,&Z

onder de voorwaarde v(x

+

pik, Z

+

qZ) I

vrna

minimali seer C I P + c24

onder de voorwaarde (3.26)

De voorwaarde VI yma is niet altijd convex (zie [CAM], [RQC]). Qm ervoor te zorgen dat dit wel convex is wordt de waarde van het operatiedoel vastgelegd,

cp+c2q =

a.

Hiermee wordt de bijbehorende voorwaarde

n n

(3 27)

welke wel convex is.

Voor een vaste waarde a kan nu worden vastgesteld of t,v< t,vmax feasible is en vervolgens kan op z'n minst een locale oplossing voor het plane search probleem

(3.26) worden gevonden.

Uit ervaringen met het werken met de verschillende methoden in de praktijk blijkt de lange-stap pad-volgende methode bij lage nauwkeurigheid langzamer te zijn dan de potentiaal reductie methoden. Bij hoge nauwkeurigheid echter is de lange-stap pad- volgende methode sneller dan de potentiaal reductie methode.

Het aantal benodigde iteraties is onafhankelijk van de grootte van het probleem. De inspanning die gelevert moet worden bij iedere iteratie is aaiankelijk van hoeveelheid structuur die in de matrix F, aanwezig is (blokdiagonaal).

3

3.5 : Phase I en Phase I-Phase 11 methoden

In het voorgaande is steeds aangenomen dat er strikt feasible primaire en duale

beginwaarden bekend waren. In de meeste gevallen is dit ook het geval, maar het komt voor dat er een of dat ze beiden niet bekend zijn. Een mogelijke oplossing voor deze situatie wordt in het onderstaande aangedragen.

De volgende gevallen zijn mogelijk :

I. Een strikt feasible x is bekend, maar geen strikt feasible 2 11. Een strikt feasible 2 is bekend, maar geen strikt feasible x

111. Een strikt feasible x noch een strikt feasible Z is bekend

Een manier om dit probleem op te lossen is met de Big-A4 methode.

Het eerste geval wordt hieronder uitgewerkt voor de andere gevallen wordt verwezen naar [BZ, 56.11.

Stel er is een strikt feasible

do)

gegeven, maar geen strikt feasible duale beginwaarde. In dit geval wordt het primair probieem (3. i) aangepast door aan het spoor vm F(x)

minimaliseer cT x

onder de voorwaarde F ( x ) 2 O

Tr F ( x ) 5 M

(3.28)

A l s Mgroot genoeg wordt gekozen zijn de oplossing van (3.28) en het origineel (3.1) hetzelfde.

Het dual probleem van (3.28) is maximaliseer

onder de voorwaarde Tr F,(Z - z l ) = c,, i = 1, ..., na - Tr F,(Z - 21) - Mz

Z2O,z20

(3.29)

Om voor dit duale probleem strikt feasible punten te berekenen is geen probleem. Eerst wordt een oplossing U = U' bepaald uit de vergelijkingen Tr <U = cz .

Neem zfo' > -rninf&In(U),O>, en stel 2fQ) = U + z@)l.

Het moeilijkste gedeelte is het bepalen van een geschikte keuze voor M.

Stel dat er optimale oplossingen x, 2, z zijn gevonden voor (3.28) en (3.29). Er geldt

TsF(x)Z

+

(A4

- Tr F(x))z = O (3.30) Er zijn nu twee mogelijkheden.

1) z = O, en er moet gelden Tr F(x) < M. In dit geval doet de extra voorwaarde in

(3.28) niet mee. Hieruit blijkt dat Mgroot genoeg is gekozen en

x

het oor- spronkelijk SP (3.1) oplost.

in dat M groter moet worden gekozen. Men kan nu met een grotere

A4

het algoritme opnieuw aanroepen.

4 : Software tj 4.1 : Inleiding

Als vervolg op de technische achtergrond van de SP methode wordt in dit hoofdstuk gekeken naar de software die gebruikt kan worden voor het oplossen van SP's.

5

4.2 behandelt het sp-pakket. Dit pakket bestaat uit de functies sp.m, bigM.m en phasel .m. Daar in de meeste gevallen zowel een primaire feasible beginwaarde als een duale feasible beginwaarde bekend is, zal alleen sp.m belicht worden. Voor de fbncties bi@d.m efi phase1.m wordt verwezen fiar [B3' $31. Vervolgens zal in

5

4.3 aandacht besteed worden aan het pakket LMitool. LMItool is een interface voor het sp-pakket. LMItool neemt een deel van het werk over dat, indien sp.m direkt aangeroepen wordt, gedaan moet worden.

§ 4.2 : Het sp-pakket

Het sp-pakket bevat software voor het oplossen van SP problemen zoals

minimaliseer C T X

onder de voorwaarde F ( x ) 4 Fo

+

x,<

+.

. .+x,F, 2 O (4.1)

De gegevens gegevens voor dit probleem zijn de vector c E

IR"

en m+l symmetrische

matrices F07..

.Jim

E R""".

Voor gebruik met Matlab zijn de functies sp.m, bigM.m en phase1.m in dit pakket aanwezig.

I. sp.m : als zowel een primaire als een duale feasible beginwaarde bekend is

11. bigM.m : als er een primaire maar geen duale feasible beginwaarde bekend is

III.

phasel .m : als zowel geen primaire als geen duale feasible beginwaarde bekend is wordt met phasel .m een primaire feasible beginwaarde bepaald. Als dit gedaan is kan verder bigM.m gebruikt worden.

54.2.1 : s p m

sp.m lost het volgende SP op

minimaliseer CT x

onder de voorwaarde F,

+

xiq

+.

. .+x,F, 2 O

en de bijbehorende duaal

maximali seer - Tr FoZ onder de voorwaarde Z = ZT 2 O

TrF,Z= ~ ~ ~ Z = l ~ . . . ~ m

(4.3)

sp.m wordt aangeroepen met :

Betekenis van de belangrijkste ingangsvariabelen

F : bevat de blokdiagonale matrices Fo, ..., F,

c : specificeerd primaire operatiedoel 80 : pïiïïìaiïe beginwaarde

ZO : duale beginwaarde

plu : bepaald de snelheid van convergentie, zie

6

3.4 voor de betekenis van v, een waarde van v tussen 10-50 wordt aangeraden

abstol : absolute convergentie

reltol : relatieve convergentie

tv: target value

maxiters : maximaal aantal iteraties

Betekenis van de belangrijkste uitgangsvariabelen

x : laatste iteratie van het primaire punt

2 : laatst iteratie van het duale punt

u1 : bevat de waarde van de primair en de duaal operatiedoel

Voor de betekenis van de andere variabelen wordt verwezen naar [B3, 53.13.

E

J 4.2.2 : Convergentie criteria

Een van de uitgangen van sp.m is een vector u1 (zie boven).

ul[O] geeft een bovengrens en ~ l [ 11 geeft een ondergrens voor de optimale waarde.

u1[0] is het primaire operatiedoel cTx voor de primaire feasible oplossing opgeslagen in

ui[ I] is het duale operatiedoel -Tr F Z voor de duale feasible oplossing opgeslagen in

z.

X

Het dualiteitsverschil is gelijk aan ul[O]-uí[l].

Er zijn vijf mogelijkheden waaronder het programma stopt :

e Het maximale aantal iteraties (maxiters) is overschreden

de absolute tolerantie (abstol) is bereikt : ul[O]-ul[l] I abstol

e de relatieve tolerantie (reltol) is bereikt :

ul[O] en ui[l] zijn positief en ul[O]-ul[l] I reltol*ul[l] of

ul[O] en ul[l] zijn allebei negatief en ul[O]-ul[1] 5 -reltol*ul[l]

de target value (tv) is bereikt : reltol is negatief en u1[0] is kleiner dan tv

e de target value (tv) is niet bereikbaar : reltol is negatief en ui[ 11 is groter dan tv

Het target value stopcriterium is erg handig voor feasibility problemen, want sp houdt

dan de gegeven tv of als er een duaal punt is waarmee het duale operatiedoel groter wordt dan tv.

3

4.2.3 : Opmerkinpen

0 De matrices F1,

...,Fm

moeten lineair onafhankelijk zijn. Zijn de matrices F,

afhankelijk dan bestaan er ;2, ongelijk aan nul waarvoor geldt

Is dit het geval dan is het probleem niet duaal feasible of kan het gereduceerd

worden tot een probleem mei mindeï vaiabelen.

0 Als sp.m tijdens de berekening wordt afgebroken is de ingangsvariabele F

verandert. F zal opnieuw moeten worden ingevoerd alvorens sp.m kan worden

aangeroepen.

0 sp.m vraagt strikt primaire en duale feasible punten. Dit betekent dat de standaard

truc om gelijkheidsvoorwaarden toe te voegen niet werkt

(Ax

= b invoeren als

m

A z q = O .

2 = 1

Ax-b 2 O, 6-Ax 2 O).

Voor informatie over de functies bigM.rn en phase1.m wordt verwezen naar [B3, $31.

5

4.3 : LMItool

LMItool is een interface voor het hiervoor behandelde sp-pakket. Zoals in het voorgaande te is zien moet voordat sp.m kan worden aangeroepen heel wat werk verricht worden. LMItool neemt dit werk voor een stuk over. LMItool is bruikbaar voor kleine en middelgrote problemen (tot een paar honderd variabelen).

LMItool lost LMI problemen (of SP problemen) op van de vorm

5

DI,.

.

.

. . ,DN

zijn de bekende matrix variabelen

5

XI,.

.

.

.

. ,Xj,,-

zijn de onbekende matrix variabelen

5

f

is een skalaire functie van Xi,..

..

. ,XN en

Dl

,...

. .

,DN , lineair in X s ,

f

is het operatiedoel

0 GI’s zijn matrix functies van XI

,...,

XH en Dl ,...,

D N

, G,’s zijn de LME’s 0

4 ’ s

zijn matrix functies van& ,..., X N en Dl ,..., DN , HJ’s zijn de LMl’s

3

4.3.1 : De functie lmitool.m

De functie lmitool genereert twee bestanden

e een solver functie welke lmisolver aanroept e een skelet van een evaluatie functie

De gebruiker dient de evaluatie functie aan te passen om het LMI probleem te definiëren waarna de solver functie kan worden aangeroepen.

Syntax : imitool(naam,varlijst,datlijst)

e naam : de naam van het op te lossen probleem

0 varlijst : bevat de onbekende matrix variabelen (gescheiden door komma’s).

e datlijst : bevat de bekende matrix variabelen (gescheiden door komma’s).

Het maximum aantal is 15. Het maximum aantal is 15.

Met het commando :

Imitool(‘ri’,’Q,Y,gam’,’A,B,C,D’)

worden door LMItool twee bestanden, ri.m (solver) en ri_eval.m (evaluatie) aangemaakt (bijlage 1).

In de solver functie zijn voor de relatieve tolerantie, absolute tolerantie, de snelheid van convergentie en het maximaal aantal iteraties reeds waarden ingevoerd. Het is aan de gebruiker om deze zo te laten of aan te passen. In de uitgewerkte voorbeelden zal

blijken dat in sommige gevallen de waarden moeten worden aangepast. .

In de evaluatie functie kan de gebruiker het LMI probleem deñniëren aan de hand van de in te voeren Lh4E’s en LMI’s. Het operatiedoel wordt gedefinieerd met OBJ.

Bij OBJ = [ ] wordt het feasibility probleem bekeken.

Bij het definieren van het probleem dient men er rekening mee te houden dat de ongelijkheden die in de evaluatie functie worden ingevoerd door LMItool gelezen worden als 2 .

Dus

A+B I O wordt geprogrammeert als -A-B 2 O

A+B > O wordt geprogrammeert als A+B-E 2 O, met O < E << 1 ,S 4.3.2 : De werking; van LMItool

De hnctie lmisolver werkt in Vier stappen

Initial p e s $ :

hier wordt voor het eerst de evaluatie functie aangeroepen, met xdatu. Uit deze eerste schatting wordt het aantal en de grootte van matrix variabelen afgeldt.

Elimination of equality constraints :

hierbij wordt herhaaldelijk de evaluatie functie aangedaan. Een kanonieke represen- tatie van de vorm

-T

minimaliseer c z

onder de voorwaarde

F,

N

+

zl&

+.

. .+z,F, 2 O,

AZ

+

b = O (4.5)

wordt geconstrueerd, waar z alle coëfficiënten van alle matrix variabelen bevat.

Elimination of variables ;

lmisolver elimineert alle overtollige variabelen. De gelijkheidsvoorwaarden worden geëlimineerd door de nulruimte N van A te berekenen en een oplossing ZO (als er een is)

van Ax

+

b = O te berekenen.

In dit stadium worden alle oplossingen van de gelijkheidsvoorwaarden gepara-

metriseerd door z =

Nx

+

zo met x een vector die de onafhankelijke variabelen bevat. Als de gelijkheidsvoorwaarden geëlimineerd zijn, wordt het probleem herformuleerd

minimaiiseer C T Z

onder de voorwaarde Fo

+

z&

+.

. .+z,F, 2 O,

met c een vector en FO,. . .,& symmetrische matrices,

x

bevat de onafhankelijke elementen van de matrix variabelen

XI,.

. . . . .Jm.

Optimization

Tenslotte roept lmisolver de functie sp.m aan om het SP probleem (4.6) op te lossen. Deze fase is onderverdeeld in een feasibility fase en een minimalisatie fase (ais er een operatiedoel gespecificeerd is). De feasibility fase wordt overgeslagen als de initial

guess feasible is.

De functie sp.m wordt aangeroepen met de optimalisatie parameters abstol, nu,

maxiters, reltol.

Bij het uitwerken van de verderop volgende voorbeelden is steeds gebruikt gemaakt

van LMztool.

In [ElG] staan enkele voorbeelden van het definiëren van verschillende regeltechnische problemen in LMI’s, zodat ze geschikt zijn om met LMItool op te lossen.

5 : Het werken met de SP methode

3

5.1 : Inleiding

In dit hoofdstuk wordt de SP methode beter bekeken aan de hand van enkele voor- beelden. In § 5.2 wordt gekeken naar het LQR probleem. Dit probleem kan zowel met klassieke methoden als met de SP methode opgelost worden. De twee methoden

worden vergeleken om zo de voor- en nadelen van de SP methode ten opzichte van de klassieke methode te bepalen. In tj 5.3 zal het inverse probleem van de optimale rege- !aar worden bekeken. at prableem is met klassieke methoden ïiet up!ûsbaaï. C e SP methode is wel geschikt om dit probleem op te lossen. De resultaten die via de SP methode verkregen zijn worden gecontroleerd met klassieke methoden.

8

5.2 : Het LOR probleem

Er wordt uitgegaan van een lineaire tijdsonaaiankelij k systeem

X = A x + B , z = C x + D u (5.1)

Het LQR probleem wordt als volgt geformuleerd : vind met beginwaarde x(0) de ingangsgrootheid u waarvoor de uitgangsenergie

1'

z'zdt geminimaliseerd wordt.

Er wordt aangenomen dat (A,&

c)

minimaal is,

d D

inverieerbaar is en d C = O.

De optimale input u kan dan worden geschreven als een toestandsterugkoppeling

u = ICY, met K = -( DTD)-' BTP.

P is de unieke positief-definite oplossing van

m

O

A'P

+

PA - PB(D~D)-'B*P

+

cTc

=

o

De optimale uitgangsenergie wordt dan gegeven door x(û)Px(û).

Met Q = P' wordt (5.2)

AQ

+

QA' - B(D~D)-'B*

+

QC'CQ =

o

_(5.3) De optimale uitgangsenergie wordt nu gegeven door x(U)Q-'x(O).

In [BI, $7.4. i] wordt gezocht naar een toestandsterugkoppelmatrix K zodat de bovengrens van de uitgangsenergie geminimaliseerd wordt. Voor lineaire tijds- onafhankelijke systemen is dit gelijk aan het minimaliseren van de uitgangsenergie. In

van K waardoor het LQR-probleem geschreven kan worden als convex probleem.

Voor een gegeven toestandterugkoppelmatrix K zal de uitgangsenergie van het sys-

teem (5. i> riet groter wordt dm x(O)Q-'x(G) wambij Q > O en Y = KQ voldoen aan

1, 97.4. i ] wordt een variabele Y geïntroduceerd. Dit leidt tot overparametrisatie

AQ

+

QAT +BY

+

YTBT (CQ

+

DY)T

Er wordt nu een toestandsterugkoppelmatrix K-YQ-’, met Y als variabele, gevonden waarbij de uitgangsenergie kleiner is dan y door het SP probleem x(O)Q-’x(û)S y en

(5.4) op te lossen. Dus minimaliseer Y onder de voorwaarde x(0)Q-’x(0) 5 y ,

L

CQ

+

DY ( 5 . 5 )

[AQ

+

QAT +BY

+

YTBT (CQ

+

DY)T

In [Bl, 57.2.31 staat dat x(û)Q-’x(O)lykan worden geschreven als

Hiermee wordt het SP probleem

minimaliseer Y

[P

;]>o

onder de voorwaarde

] < O

AQ

+

QAT +BY

+

YTBT (CQ

+

DY)T

CQ +DY -I

met Q > O en Y =KQ.

Met LMItool wordt nu een oplossing voor dit probleem gezocht. De werking van LMItool is reeds in het voorgaand hoofdstuk behandelt. De voor dit probleem gebruikte LPJE’S en LMI’s zijn

1 2 0

AQ

+

QAT +BY

+

YTBT (CQ +DY)T

CQ

+

DY -I

De eis Q = Q’ moet hier extra toegevoegd worden. In de theorie volgt deze eis

automatisch uit Q > O (zie Notaties).

ydient geminimaliseerd te worden dus is yde objective. De bij dit probleem behorende

Voorbeeld 5.1 :

Bij het eerste voorbeeld is het onderstaande systeem gekozen (bijlage 3) waarbij de matrices C en D moeten voldoen aan DTC = O en

DTD

moet inverteerbaar zijn.

( 5 . 8 )

-

2 1

Er wordt uitgegaan van de situatie dat geen enkele indicatie is van wat de uitkomst van

Q en Y zal zijn. Daarom is als beginvoorwaarde voor Q een nulmatrix en voor Y een nulvector gekozen. Als beginvoorwaarde voor het operatiedoel yis gekozen yb = lo4. Nu woïdt de solver hnctie (zie bijlage 2) aangeroepen. De gevonden waarden VOOT Q,

Y en y staan in bijlage 4. Met P = Q-' en K = Ye-' volgt

53 1339 1339 6149.1

K = [-1.1272 -51.76031 ; P = _(5.9)

Met Matlab wordt hetzelfde LQR probleem opgelost met een bestaande methode (NB : matlab gebruikt voor het LQR probleem u = -Kx terwijl in dit verhaal wordt

uitgegaan van u=Kx de K die matlab vindt dient met min te worden vermenigvuldigd voor het wordt gebruikt om het bovenstaande LMI probleem op te lossen)

Dit levert K* en P* (zie bijlage 4)

(MY).

K* = [-1.1272 -51.76031 ; P* =

[

53

1339 6149.1 (5.10)

Het verschil tussen K, K* is van de orde lom5 en die tussen P, P* is van de orde

Hieruit is te concluderen dat beide methode niet voor elkaar onderdoen wat betreft het resultaat.

Echter als de numerieke nauwkeurigheid bij de LMI[ methode te klein wordt gekozen stopt het programma de berekening en worden er geen resultaten verkregen.

Deze numerieke nauwkeurigheid wordt in de solver functie (zie bijlage 2 ) opgegeven

middels de relatieve en de absolute tolerantie. In het bovenstaand voorbeeld is voor

deze tolermties een w a r d e genomen wordt

genomen wordt het programma in de optimaliseringsfase onderbroken met als fout- melding dat er niet aan de feasibility eisen is voldaan.

.

Dit probleem treedt ook op wanneer er bij eenzelfde tolerantie de matrix D vergroot wordt. Bij vergroting van de nauwkeurigheid zal het programma dan echter wel een oplossing geven.

Voorbeeld 5.2 :

Om het verschil in oplostijd tussen beide methoden te vergelijken wordt nu een groter systeem (zie bijlage 5 ) bekeken

X=P,Y+Ru ; z = C x + D u O 1 0 O O O -9.81 -0.8918 O O -0.0909 9.0909

c =

Jz

; D = (5.11)

Ook hier zijn C en D weer zo gekozen dat aan de eisen D'C = O, 0% inverteerbaar

voldaan is.

Net als bij het voorgaande voorbeeld wordt voor de beginschatting van Q een

nulmatrix en voor die van Y een nulvector genomen. De beginschatting voor y is 1 O4 . Nu wordt de solver functie (zie bijlage 2) aangeroepen. De gevonden waarden van Q,

Y en yzijn te vinden in bijlage 6. Met P = Q-' e n K = Ye'' volgt hieruit

K = [1.1029 1.7299 -8.7322 -1.69911

r

5.4764 1.9926 -1.4366 -0.2426 I 1,9923 2.7293 -7.6707 -0,3805 - 1.4366 -7.6707 28.1 1 54 1.92 1 O -0.2426 -0.3805 1.9210 0.3738 (5.12)

Wederom worden deze waarden vergeleken met de waarden die een reeds bestaande oplosmethode geeft (bijlage 6 ) . Voor K* en P* volgen hieruit

K* = [I1029 1.7299 -8.7321 -1.69911 5.4764 1.9926 -1.4366 -0.2426 1.9926 2.7293 -7.6707 -0.3 805 -1.4366 -7.6707 28.1 154 1.9210 -0.2426 -0.3805 1.9210 0.3738 (5.13)

De verschillen tussen K, K* en tussen P, P* zijn ook in dit voorbeeld klein dus wat dat aangaat doen beide methoden niet veel voor elkaar onder. Het verschil tussen de methoden zit ‘m meer in de benodigde opiostijd. -Waar de idassieke methode (matlab functie lqry) in een oogwenk het resultaat geeft, heeft de LMI methode aanmerkelijk meer tijd nodig (-20 seconden) en dan hebben nog pas over een systeem met dimensie 4. Bij grotere systemen lopen de verschillen nog verder op.

,a

5.3 : Het inverse probleem van de optimale regelaar

Normaal wordt bij een optimale regelaar aan de hand van weegmatrices Q en I! een toestandstenigkoppelmatrix K bepaald, zodat de regelkring (figuur 3) gestabiliseerd wordt.

toestandsvergelijking

i

= A x + B u

I-

optimale lineaire regeIwet

u

=

Kx

figuur 3

Dit probleem wordt als volgt geformuleerd [B 1,

0

10.61: Gegeven een systeem

X = A X + B î .

r i

met x(O) =XO , (A,B) stabiliseerbaar, @,A) detecteerbaar en R > O. Kan de oplossing

van het LQR probleem (zie

5

5.2) geschreven worden als een toestandsterugkoppeling

u=Kx met K=-R-’BTP, waarbij P de unieke niet-negatieve oplossing is van

ATP

+

PA - PBR-’BTP

+

Q = O (5.15)

met weegmatrices Q en R.

Bij het inverse probleem van de optimale regelaar (de naam zegt het al) wordt het pro-

bleem omgedraaid. Ook hier wordt m~ een temgkoppehatnk K m weegmatrices Q

en R gezocht zodat de regelkring (figuur 3) gestabiliseerd wordt.

Neem aan dat de terugkoppelmatrix K bekend is (bv. via poolplaatsing). Er wordt nu gezocht naar geschikte weegmatrices Q 2 O en R > O, zodat (Q,A) detecteerbaar is en u =

lix

de optimale regelaar voor het bijbehorende LQR probleem is.

De oplossing van het inverse probleem wordt gevonden door te zoeken naar een Q 2 O

en R > O, zodat er een niet-negative P en een positief-definite PI bestaan die voldoen aan

( A

+

B K ) ~

P

+

P ( A

+

B K )

+

K ~ R I ~ T

+

Q = O, B T P + K =

o

(5.16) en A T P , + < A < 0

De voorwaarde met PI is hetzelfde als de eis dat @,A) detecteerbaar moet zijn. Met behulp van LMItool wordt nu een oplossing gezocht voor dit probleem. De werking van LMItool is reeds in een voorgaand hoofdstuk behandelt. Hieronder wordt aangegeven welke LME’s en LMI’s gebruikt worden.

LMEi’s : ( A

+

BTP

+

RK

=

o,

P

+

P(A

+

B K )

+

KTRK

+

Q =

o,

P = P T , Q = Q T e n < = < T (5.17)

Voor de laatste drie voorwaarden geldt hetzelfde als voor de voorwaarde QT= Q bij het LQR probleem (zie

0

5.2).

LMI’S :

-ATP,-&A+Q>O,

Q 2 0 , R>0, P > O e n < > O (5.18)

Daar we nu een feasibility probleem bekijken is er geen operatiedoel (OBJ=[])

De bij dit probleem behorende solver functie en de aangepaste evahatie functie zijn te

Voorbeeld 5.3 :

Bij dit voorbeeld wordt het volgende systeem bekeken (zie bijlage 5 )

i ( t ) = Ax(t)

+

h ( t ) A = O 1 0 O O O -9.81 O O O 0 1 -8.89’18

o

0 -0.0909 ; B = (5.19)

Via matlab (lqr-methode, zie opmerking bij voorbeeld 5.1) wordt nu een terugkoppel- matrix K bepaald (normalitair is het de bedoeling dat K gevonden wordt door middel van poolplaatsing) door het LQR probleem op te lossen met voor Q e en Re

[l

o o

o1

Q e = o

/ O

o o

1 O I ; R e = i

o

(5.20)

Met deze waarden volgt voor de terugkoppelmatrix K

K = [1.1029 1.5492 -6.8349 -1.57231 (5.21)

Door nu met deze terugkoppelmatrix K het LMí probleem op te lossen zouden er een

Q en R moeten volgen die hetzelfde systeem beschrijven.

Dit is te controIeren door met deze, via LMítool berekende, Q en R het LQR probleem op te lossen, er moet dan eenzelfde toestandsterugkoppelmatrix

R

volgen.

Controle van de uitkomst is niet mogelijk door Q, R te vergelijken met Q e , Re.

Q en R hoeven namelijk niet gelijk te zijn aan Qe en Re om toch hetzelfde systeem te beschrijven.

Met deze A, B, K, nulmatrices als beginvoorwaarden voor Q,

P,

PI

en nul als beginvoorwaarde voor R wordt nu de solver functie (zie bijlage 7) aangeroepen. De gevonden waarden voor de matrices

P,

PI,

Q en R staan in bijlage 8.

Om dit resultaat te controleren wordt met deze waarden voor Q en R opnieuw het LQR probleem opgelost om een terugkoppelmatrix

If

te bepalen.

K* = [l,1029 1,5492 -6,8349 -1,57231 (5.22)

Als we de waarde van K* vergelijken met de waarde van K blijkt er een verschil tussen te zitten ter grootte van Hieruit blijkt dat bij het oplossen van het LMI probleem, beschreven in 5.17 en 5.18, er een Q en R volgen die voldoen aan de gestelde eisen. De SP netkmde is dus geschikt VOW het oplosser, vm het iriverse probleem vm de

Voorbeeld 5.4 :

Bij dit voorbeeld wordt gekeken naar het volgende systeem (bijlage 3)

i ( t ) = h ( t )

+

Bzt(t)

Via matlab wordt een matrix K bepaald door het LQR probleem op te lossen . Voor Q, en R, zijn de volgende waarden genomen

1 0

Q,

=[o

Re = 1

Met deze waarden volgde voor de terugkoppelmatrix K K =

[-

1,127 1 -5 1,76031

(5.23)

(5.24)

(5.25) Er geldt &er natuurlijk hetzelfde als voor voorbeeld 5.3. Door nu met deze

terugkoppelmatrix K het LMI probleem op te lossen zouden er een Q en R moeten volgen die hetzelfde systeem beschrijven. Met deze Q en R moet bij het oplossen van het LQR probleem dezelfde terugkoppelmatrix K moeten volgen.

Voor de beginvoorwaarden van P, Pr en Q zijn wederom nulmatrices gekozen en voor de beginvoorwaarde van R is nul gekozen. Met deze waarden wordt de solver functie aangeroepen (zie bijlage 7).

De gevonden waarden voor

P,

P I , Q en R zijn te vinden in bijlage 9.

Door deze berekende Q en R wordt nu op dezelfde manier als bij voorbeeld 5.3 een terugkoppelmatrix K* gevonden

K* = [-1.1688 -57.38571 (5.26)

Bij het vergelijken van K en K* blijkt dat er een behoorlijk verschil bestaat tussen de terugkoppelmatrix

x"

, welke is verkregen door het LQR probleem op te lossen gegeven de met de SP methode berekende Q en R, en K.

Ter controle is gekeken of in dit geval voor de berekende P, PI, Q en R voldaan is aan de bij het inverse probleem gestelde eisen. In bijlage 10 zijn de resultaten van deze controle opgenomen. Uit deze resultaten blijkt dat aan alle gestelde voorwaarden wordt voldaan.

Het verschil tussen K en K* kan ook ontstaan zijn doordat de numerieke

nauwkeurigheid van het algoritme niet toereikend is. De relatieve en de absolute tolerantie worden gedefinieerd in de solver functie (zie bijlage 11, relltol en abstol).

Ter controle zijn beide verkleint met een factor 1000.

K* = [-1.1688 -57.38551 (5.27)

Deze aanpassing sorteert niet het gewenste succes. Voor dit voorbeeld lijkt de SP methode (in deze vorm) niet te voldoen.

Voorbeeld 5.5 :

Dit voorbeeld is bijna hetzelfde als bovenstaand voorbeeld. Alleen ingangsmatrix B is een factor 1

o4

groter genomen

B = [ 8.4 O 175

1

(5.28)

verder gelden dezelfde waarden als bij het voorgaande voorbeeld.

Met de matrix B volgt voor de oplossing van het LQR probleem met Qe , Re uit (5.24)

K = 1-1.0001 -1.83741 (5.29)

De SP methode geeft, met dezelfde beginwaarden als bij voorbeeld 5.4, een Q, R (zie

bijlage 13) waarbij voor de oplossing van het LQR probleem, met de bestaande methode, volgt

K* = [-1.0001 -1.8374] (5.30)

Het verschil tussen K en K* is ter grootte dit voorbeeld wel voldoet.

Het enige verschil tussen voorbeeld 5.4 en voorbeeld 5.5 is dat de ingangsmatrix B een factor lo4 groter is. Waarom met de SP methode bij voorbeeld 5.5, in tegenstelling tot voorbeeld 5.4, wel een goede oplossing gevonden is, is niet goed verklaarbaar.