EINDWERK: Actieve demping van een ski.
Studiegebied Industriële Wetenschappen en Technologie Opleiding Elektromechanica Optie Elektromechanica Academiejaar 2005-2006
Karel Deprez
EINDWERK: Actieve demping van een ski.
Studiegebied Industriële Wetenschappen en Technologie Opleiding Elektromechanica Optie Elektromechanica Academiejaar 2005-2006
Karel Deprez
Voorwoord.
Om te beginnen wil ik iedereen bedanken die de realisatie van dit eindwerk mogelijk heeft gemaakt.
In de eerste plaats wil ik Dhr. dr. ir. Frederik D’hulster bedanken voor de begeleiding van dit eindwerk. Zijn vele tips en aanwijzingen hebben een zeer positieve invloed gehad op dit eindwerk. Ook Dhr. dr. ing. Kurt Stockman wil ik bedanken voor de begeleiding tijdens het maken van de regelaar.
Het departement PIH van de Hogeschool West-Vlaanderen bedank ik voor de middelen die me werden aangeboden om mijn eindwerk te kunnen uitvoeren.
Tot slot wil ik mijn dank betuigen aan mijn ouders voor hun altijd aanwezige steun.
Karel Deprez.
I
Inhoudsopgave. inleiding ...............................................................................................................................1 1 Hoofdstuk 1: Doelstellingen. .......................................................................................3 2 Hoofdstuk 2: Proefopstelling. ......................................................................................4 2.1 Totale opstelling...................................................................................................4 2.2 Skilat ....................................................................................................................5 2.3 Actuator................................................................................................................5 2.4 Versnellingsopnemer ...........................................................................................6 2.5 Meetversterker en power amplifier......................................................................6 2.6 De BNC box.........................................................................................................7 2.7 dSPACE Controldesk®. ......................................................................................7 2.8 Matlab/Simulink®. ..............................................................................................7 3 Hoofdstuk 3: Trillingen................................................................................................8 3.1 Problematiek van trillingen..................................................................................8 3.2 Transversale trilling van een balk aan één zijde ingeklemd. ...............................8 4 Hoofdstuk 4: Voorstellen van systemen[7]................................................................12 4.1 Voorstelling in het s-domein..............................................................................12 4.1.1 Laplace transformatie.................................................................................12 4.1.2 Betekenis van nullen en polen ...................................................................13 4.1.3 Voorstellen in het frequentiedomein..........................................................14 4.2 Voorstelling in het Z-domein.............................................................................16 4.2.1 Algemeen ...................................................................................................16 4.2.2 De Z – transformatie ..................................................................................17 4.2.3 Betekenis van nullen en polen ...................................................................18 4.2.4 Blokschema van een digitaal regelschema. ...............................................18 5 Hoofdstuk 5: Systeem identificatie, theoretische benadering.[8] ..............................20 5.1 Algemene theorie. ..............................................................................................20 5.1.1 De signalen.................................................................................................20 5.1.2 Polynoomvoorstelling van een transferfunctie. .........................................21 5.1.3 Toestandsruimtevoorstelling van een transferfunctie. ...............................22 5.2 Modellen voor het identificatieproces. .............................................................22 5.2.1 ARX ...........................................................................................................22 5.2.2 N4SID ........................................................................................................24 5.3 Graphical User Interface ....................................................................................24 5.4 Datavoorstelling.................................................................................................25 5.4.1 Data importeren .........................................................................................26 5.4.2 Data bestuderen..........................................................................................27 5.4.3 Voorbewerking van data ............................................................................28 5.4.4 Modellen schatten m.b.v. data ...................................................................28 5.4.5 Modellen analyseren ..................................................................................30 6 Hoofdstuk 6: Identificatie van een 2de orde systeem.[6]............................................31 6.1 Meetopstelling....................................................................................................31 6.2 Instelling. ...........................................................................................................32 6.3 Identificatie. .......................................................................................................32 6.3.1 ARX ...........................................................................................................33 6.3.2 N4sid ..........................................................................................................36 7 Hoofdstuk 7: Identificatie van een 1e orde systeem...................................................39 7.1 Instelling ............................................................................................................39
II
7.2 Identificatie. .......................................................................................................40 7.2.1 ARX ...........................................................................................................41 7.2.2 N4SID ........................................................................................................43 7.2.3 Identificatie d.m.v. GUI .............................................................................45 8 Hoofdstuk 8: Identificatie van de skilat .....................................................................48 8.1 Systeembeschrijving. .........................................................................................48 8.2 Identificatie ........................................................................................................49 8.2.1 Identificatie m.b.v. de GUI. .......................................................................50 9 Hoofdstuk 9: Analyse van een toestandsruimte.[9] ...................................................55 9.1 Algemene voorstelling .......................................................................................55 9.2 Wiskundige modellen van lineaire systemen.[1] ...............................................56 9.3 Controleerbaarheid[2] ........................................................................................57 9.4 Observeerbaarheid[2].........................................................................................58 9.5 Stabiliteit.[7] ......................................................................................................58 9.6 Poolplaatsing door toestandsterugkoppeling. ....................................................59 9.7 Toestandschatter. ...............................................................................................63 10 Hoofstuk 10: Ontwerp toestandsruimteregelaar.[7]...................................................65 10.1 Toestandsterugkoppeling. ..................................................................................67 10.2 Toestandsschatter...............................................................................................71 10.3 Toestandsruimteregelaar[5] ...............................................................................74 11 Hoofdstuk 11: Actieve demping van de ski...............................................................77 11.1 Implementatie van de regelaar. ..........................................................................77 11.2 Controle werking. ..............................................................................................78 besluit.................................................................................................................................80
III
Lijst symbolen en afkortingen. f: F: B: I: L: M: V: E: I: ρ: A: a: G(z): H(z): Y(s) R(s) K ω: φ: Ts: δ (t ) : u(t): y(t): e(t): D: τ: x u y λ: I: L: TF: GUI:
frequentie [Hz] kracht [N] inductie [T] stroom [I] lengte [m] moment [Nm] schuifkracht [N] elasticiteitsmodulus[N/m²] traagheidsmoment [m4] massadichtheid[kg/m³] oppervlakte[m²] versnelling [m/s²] tranferfunctie van het systeem verband tussen uitgang en storing uitgang analoog systeem ingang analoog systeem versterkingsfactor hoekfrequentie [rad/s] fasehoek [°] sampletijd dirac impuls ingang uitgang ruis demping tijdsconstante de toestandsvector de ingangsvector de uitgangsvector de karakteristieke wortels of eigenwaarden eenheidsvector terugkoppelmatrix transferfunctie graphical user interface
IV
Figurenlijst. fig. 2.1: Totale proefopstelling van de skilat. ..................................................................... 4 fig. 2.2: De actuator en sensor van de opstelling. ............................................................... 6 fig. 3.1: Eénzijdig ingeklemde balk. ................................................................................... 8 fig. 3.2: Segmentje van de éénzijdig ingeklemde balk. ...................................................... 9 fig. 4.1: Algemene voorstelling van een transferfunctie van een proces. ......................... 12 fig. 4.2: Nulpunten-polen-diagram van een 1e orde systeem. ........................................... 13 fig. 4.3: Verband tussen de ligging van de pool en het transiënt gedrag. ......................... 14 fig. 4.4: Algemene voorstelling van een bodekarakteristiek............................................. 15 fig. 4.5: Relatie tussen tijdsgrafiek, nulpunten-polen-diagram en de bodediagram. ........ 15 fig. 4.6: AM en FM aflezen uit een bodekarakteristiek. ................................................... 16 fig. 4.7: Overgang van analoog naar digitaal signaal........................................................ 17 fig. 4.8: Relatie tussen het Z-domein en het S-domein..................................................... 18 fig. 4.9: Blokschema van een digitale regelaar. ................................................................ 19 fig. 5.1: Algemene signaalvoorstelling met u=input,y=output en e=ruis. ........................ 20 fig. 5.2: Algemene polynoomvoorstelling ........................................................................ 21 fig. 5.3: ARX voorstelling. ............................................................................................... 23 fig. 5.4: Graphical user interface die de interactie met de toolbox verzorgt..................... 25 fig. 5.5: Data importeren in de GUI.................................................................................. 27 fig. 5.6: Data bestuderen in de GUI. ................................................................................. 28 fig. 5.7: Gewenst model met zijn orders ingeven. ............................................................ 29 fig. 5.8: Gelijkenis controleren tussen model en werkelijke signaal. ............................... 30 fig. 6.1: Meetopstelling voor de identificatie van een gekend systeem. ........................... 31 fig. 6.2: Bodekarakteristiek van het 2de orde systeem....................................................... 33 fig. 6.3: Voorstellen van ingang en uitgang van het 2de orde systeem.............................. 34 fig. 6.4: Gelijkenis controleren tussen model en werkelijk signaal. ................................. 34 fig. 6.5: Bodeplot van het 2de orde model. ........................................................................ 36 fig. 6.6: Gelijkenis controleren tussen model en werkelijk signaal. ................................. 37 fig. 6.7: Bodeplot van 2de orde model. .............................................................................. 38 fig. 7.1: Ingang en uitgang, voorstelling van 1e orde systeem. ......................................... 40 fig. 7.2: Bodeplot van het 1e orde systeem. ...................................................................... 41 fig. 7.3: Gelijkenis controleren tussen model en werkelijk signaal. ................................. 42 fig. 7.4: Bodeplot van 1e orde model. ............................................................................... 43 fig. 7.5: Gelijkenis controleren tussen model en werkelijk signaal. ................................. 44 fig. 7.6: Data importeren via deGUI. ................................................................................ 45 fig. 7.7: Parametric models,keuze van het model. ............................................................ 45 fig. 7.8: Gelijkenis controleren tussen model en werkelijk signaal. ................................. 46 fig. 7.9: Info opvragen van gewenste model..................................................................... 46 fig. 7.10: Eénheidscirkel van het gekozen model. ............................................................ 47 fig. 8.1: Proefopstelling met actuator en sensor................................................................ 48 fig. 8.2: In- en uitgang van aangelegd sweep-signaal....................................................... 49 fig. 8.3: In- en uitgang van aangelegd sweep-signaal in de GUI...................................... 50 fig. 8.4: Parametric models,keuze van gewenst model..................................................... 51 fig. 8.5: Gelijkenis controleren tussen model en werkelijk signaal. ................................. 51 fig. 8.6: Voorstelling van de éénheidscirkel van het model.............................................. 52 fig. 8.7: Eénheidscirkel model met dominante polen en nullen (ingezoomd). ................. 53 fig. 8.8: Bodeplot van de skilat bekomen met sweep-signalen......................................... 54 fig. 9.1: Blokschema toestandsvergelijking. ..................................................................... 56 fig. 9.2: Toestandsterugkoppeling..................................................................................... 60
V
fig. 9.3: Vlugge demping met behulp van het nulpunten-polen-diagram. ........................ 61 fig. 9.4: Algemene voorstelling van toestandsregelaar met toestandsschatting................ 63 fig. 10.1: Controle toestandsmodel. .................................................................................. 66 fig. 10.2: Scoopbeeld van controle toestandsmodel.......................................................... 67 fig. 10.3: Nulpunten-polen-diagram van de skilat. ........................................................... 68 fig. 10.4: Schema toestandsterugkoppeling. ..................................................................... 68 fig. 10.5: Scoop 1: Voorstelling van een gedempt systeem.............................................. 69 fig. 10.6: Scoop 1: Voorstelling een snel gedempt systeem. ............................................ 70 fig. 10.7: Scoop 1: Voorstelling van een zeer snel en gedempt systeem. ......................... 71 fig. 10.8: Toestandschatter................................................................................................ 72 fig. 10.9: Scoop1: Uitgangssignaal van de toestandsschatter. .......................................... 73 fig. 10.10: Totale toestandsruimteregelaar........................................................................ 74 fig. 10.11: Scoop1: Demping van signaal......................................................................... 75 fig. 10.12: Scoop2: Stuursignaal van de regelaar. ............................................................ 76 fig. 11.1: Toestandsruimteregelaar van de skilat. ............................................................. 77 fig. 11.2: Algemene voorstelling van de regelaar ............................................................. 78 fig. 11.3: Vergelijking gedempte en niet-gedempte trilling.............................................. 78
VI
INLEIDING Trillingen zijn in de meeste gevallen storend. Denken we maar aan het nauwkeurig positioneren waar ongewenste trillingen voor grote problemen kunnen zorgen. Kleine trillingen kunnen soms andere onderdelen aanstoten en zo toch grote gevolgen hebben. Daarom is het van uiterst groot belang om trillingen, hoe klein ook, te kunnen opmeten, analyseren en compenseren.
Ook bij het skiën kunnen trillingen in de skilat zeer storend zijn. Hierdoor zal het contact met de sneeuw verminderen. Dit zal negatieve gevolgen hebben voor de bestuurbaarheid van de skilat waardoor de snelheid lager zal liggen.
Trillingen kunnen op twee manieren tegengegaan worden. Een éérste manier wordt de passieve trillingsonderdrukking genoemd. Deze bestaat erin om de eigenschappen van de structuur van het materiaal te wijzigen. Dit kan door bijvoorbeeld de massa of de stijfheid te veranderen. Het is niet altijd mogelijk om dit toe te passen. De tweede methode wordt de actieve trillingsonderdrukking genoemd. Deze methode zal met actuatoren en sensoren de trilling controleren en onderdrukken. Aan de structuur van de skilat kan niets meer veranderd worden. Daarom wordt deze laatste methode in dit eindwerk toegepast.
Het doel van dit eindwerk is het dempen van trillingen bij een skilat. Daarvoor moet eerst een identificatie van de skilat plaatsvinden. Met behulp van een nauwkeurig wiskundig model van deze skilat is het mogelijk om een regelaar te ontwerpen die de trillingen actief dempt.
Dit eindwerk is in 11 hoofdstukken onderverdeeld.Om te beginnen wordt de problematiek van trillingen in een skilat besproken. Het tweede hoofdstuk bespreekt de testopstelling en de gebruikte software. Het derde hoofdstuk schetst een algemeen theoretisch beeld van trillingen. In hoofdstuk 4 wordt een korte voorstelling gegeven van de verschillende systemen terwijl hoofdstuk 5 een uitgebreide theoretisch studie geeft omtrent identificatie. In hoofdstukken 6 en 7 wordt er een identificatie toegepast op een gekend systeem zodat in hoofdstuk 8 dan vlot de ski kan geïdentificeerd worden. Hoofdstuk 9 vertelt alles omtrent de theorie van toestandsruimteregelaars zodat in
1
hoofdstuk 10 dit praktisch kan uitgewerkt worden. In het laatste hoofdstuk wordt dan uiteindelijk de regelaar geïmplementeerd in het systeem en gevalideerd.
2
1 HOOFDSTUK 1: DOELSTELLINGEN. Mensen komen in hun vrije tijd steeds meer in beweging en zoeken de spanning op. Dit uit zich onder andere in het beoefenen van dynamische sporten zoals skiën. Bij deze sport wordt het uiterste van het materiaal gevergd. Veelal zullen afdalingen op harde sneeuw en ijs trillingen veroorzaken. Kleinere ski’s en harde sneeuwoppervlakten zorgen voor een nog groter trillingsprobleem. Door het trillen van de skilat vermindert het contact van de zijkant van de ski met de sneeuw.Bij het skiën geldt het volgende: geen contact betekent geen controle. Door het verminderen van het contact tussen ski en sneeuw worden de besturingsmogelijkheden van de skiër beperkt. Dit kan zowel gevaarlijk zijn voor de skiër als voor zijn onmiddellijke omgeving. Het tegengaan van trillingen in de skilat gaat gepaard met enkele moeilijkheden. De trillingen in de skilat zullen constant veranderen tijdens de afdaling. Het parcours, de beweging en de oppervlaktegesteldheid van de sneeuw zijn nooit constant. Daarom is het bijna onmogelijk om een trillingscontrole te creëren om een bepaald type van trilling tegen te gaan. De trillingen zullen immers onvoorspelbaar veranderen tijdens iedere nieuwe afdaling van de berg. Om deze trillingen beter te onderzoeken wordt er een proefopstelling opgebouwd. Het is van zeer groot belang om te weten hoe de ski zal reageren op veranderende trillingen. Daarom zal de skilat geïdentificeerd moeten worden. Er kan pas aan demping worden gedacht, éénmaal de identificatie voltooid is. Actieve demping komt steeds meer voor bij skilatten. Steeds vertoont de actieve demping volgend schema: een sensor die de trillingen zal opmeten en een actuator die deze trillingen zal tegen gaan. Ook bij de opgebouwde opstelling wordt dit principe toegepast.
3
2 HOOFDSTUK 2: PROEFOPSTELLING. Het doel van dit eindwerk is het wegwerken van ongewenste trillingen op een skilat. Daarvoor moet een identificatie van de skilat gebeuren. De trillingen worden weggewerkt door middel van een toestandsruimteregelaar. Deze regelaar meet met behulp van een versnellingsopnemer de ongewenste trilling op en zal vervolgens een actuator aansturen die deze trilling dan meteen gaat tegenwerken. Dit hoofdstuk beschrijft de proefopstelling alsook de gebruikte software voor het binnenlezen en verwerken van de meetsignalen.
2.1
Totale opstelling.
Figuur 2.1 toont de proefopstelling van dit eindwerk.
fig. 2.1: Totale proefopstelling van de skilat. 1) De skilat
4) De meetversterker
2) De actuator (trillingsopwekker)
5) Power amplifier
3) De sensor (versnellingsopnemer)
6) BNC-box
4
2.2
Skilat
Vele ski’s hebben een kern van hout. Deze kern wordt dan bedekt door een versterkende koolstofvezel. Er worden nu ook ski’s gebouwd die volledig van koolstofvezels gemaakt zijn. De randen worden dan meestal met een metalen legering afgewerkt. De eigenfrequenties worden via metingen achterhaald. Het analytisch of het numeriek bepalen via een eindige elementen software is niet mogelijk. De oorzaak hiervan is dat de skilat geen rechte vorm heeft en dus zeer moeilijk te tekenen is. Dat de skilat uit meerdere lagen van verschillende materialen opgebouwd is, bemoeilijkt het natuurlijk ook. Bij de meting wordt de uitwijking bekeken voor een frequentie tot 50 Hz. De 2 eigenfrequenties die gevonden zijn bij de meting: f1= 12.5 Hz f2= 42.5 Hz Dit is terug te vinden in de bijlage pagina 2.
2.3
Actuator.
De actuator wordt op de grond geplaatst en de triltafel wordt vastgemaakt aan de skilat. Als de actuator een trilling opwekt, dan zal er een kracht op de skilat worden uitgeoefend. Doordat hij op de grond staat, zal hij dan een even grote reactiekracht op dit oppervlak uitoefenen. De actuator heeft binnenin een permanente magneet met daarrond een spoel, die zich dus steeds in het magnetisch veld ervan bevindt. Wanneer er dus een stroom door de spoel wordt gestuurd, wekt dit een kracht op volgens: F=BIL
(2.1)
Met B: de magnetische inductie.[T] L: de lengte van de stroomvoerende geleider in het magnetisch veld.[m] I: de stroom door de spoel.[A]
In de opstelling wordt volgende actuator gebruikt: (Figuur 2.2 toont de actuator) Merk: Gearing & Watson Type : V40 Gegevens zijn terug te vinden in de bijlage pagina 3.
5
fig. 2.2: De actuator en sensor van de opstelling.
2.4
Versnellingsopnemer
De versnellingsopnemer die gebruikt wordt is een piëzo-elektrisch type. Dit type wordt momenteel het meest gebruikt voor trillingsmetingen. Het hart van een piëzo-elektrische versnellingsopnemer is een schijfje piëzo-elektrisch materiaal, meestal een kunstmatig gepolariseerd keramisch materiaal. Wanneer het mechanisch wordt belast, hetzij in druk, trek of schuifrichting, verschijnt er een elektrische lading op de vlakken van het kristal, welke evenredig is met de inwerkende kracht. Wanneer er op de opnemer een trilling inwerkt, dan oefent de massa een kracht uit op het piëzo-elektrisch element. Deze is evenredig met de versnelling want F = m . a . Deze versnellingsopnemer wordt meestal in 2 uitvoeringen gebruikt: het drukbelaste type, waarmee de massa een drukkracht uitoefent en een schuifbelast type waarbij de massa een schuifkracht uitoefent. Gegevens zijn terug te vinden in de bijlage pagina 4 en 5.
2.5
Meetversterker en power amplifier.
De meetversterker verzorgt de voeding voor de versnellingsopnemer. Daarbij zal ook het signaal afkomstig ervan versterkt worden. Bij de versterker kan er 10x versterkt worden. Dit is voor de opstelling zeker en vast nodig. Dit komt door het feit dat het versnellingssignaal zeer klein is. Dit signaal wordt dan via de dSPACE-kaart binnen gelezen. Bij dSPACE wordt het analoog signaal omgezet naar bits, maar omdat het signaal te klein was zonder versterking was het verkregen signaal zeer hoekig en zeker en
6
vast niet mooi sinusoidaal. dSPACE werkt met een signaal tussen 0 en 1, daarvoor heeft het een bepaald aantal bits ter beschikking. Doordat het signaal zonder versterking zeer klein was, werd er maar een klein aantal bits gebruikt om het signaal weer te geven. Door de versterking wordt het aantal gebruikte bits dus meteen de hoogte in gedreven. Dat heeft tot resultaat een mooier signaal, waarmee dan makkelijker kan verder gewerkt worden. Ook wordt er gebruik gemaakt van de ruisfunctie, die toelaat om ruis te filteren. De power amplifier verstuurt het gepaste signaal naar de actuator.
2.6
De BNC box
Om op een eenvoudige manier verbindingen te kunnen maken tussen de ingangen en uitgangen van de dSPACE DS1104-kaart en de fysische signalen wordt er gebruik gemaakt van een BNC box. Er kan een analoog signaal digitaal binnengenomen of uitgestuurd worden met dSPACE en omgekeerd.
2.7
dSPACE Controldesk®.
Het programma dSPACE ControlDesk® zorgt voor een user-interface en visualiseert de aangelegde trilling(sinus) en het outputsignaal a(t). Zo kunnen de opgemeten signalen gevisualiseerd en verwerkt worden. De dSPACE-kaart in een PCI slot van de PC heeft haar eigen processor zodat de processor van de PC ontlast wordt en zo een constante sampletijd gegarandeerd is. Met hulp van Matlab simulink® wordt een model gemaakt dat van toepassing is op de opstelling. Bij de opstelling wordt gebruik gemaakt van een BNC box zodat op die manier verbindingen kunnen gemaakt worden tussen de ingangen en uitgangen van de dSPACE-kaart en de signalen.
2.8
Matlab/Simulink®.
Matlab/Simulink® is een software pakket dat de gebruiker toelaat op een relatief eenvoudige manier diverse numerieke berekeningen uit te voeren. Daarnaast heeft Matlab/Simulink® zeer veel mogelijkheden voor visualisatie van functies en data, wat kan bijdragen tot een beter begrip. Voor dit eindwerk was Matlab/Simulink® zeer belangrijk, het is zowel gebruikt voor de identificatie van de signalen als voor het maken van de regelaar. Dit wordt in volgende hoofdstukken nog uitgebreid besproken.
7
3 HOOFDSTUK 3: TRILLINGEN. 3.1
Problematiek van trillingen.
Bij snelle afdalingen op harde sneeuw en ijs gaan ski’s trillen. Door de trilling vermindert het contact tussen de ski en de sneeuw zodat de besturing ervan moeilijker wordt. Dit heeft een negatief effect op de beheersbaarheid ervan waardoor de snelheid wordt beperkt. Om een beter beeld te krijgen van de trillingen in de ski wordt de transversale trilling van een balk, aan één zijde ingeklemd, besproken. De inklemming kan overeenkomen met de bevestiging van de ski ter hoogte van de schoen.
3.2
Transversale trilling van een balk aan één zijde ingeklemd.
De dynamica van een éénzijdig ingeklemde balk wordt hier kort besproken.
fig. 3.1: Eénzijdig ingeklemde balk.
Er wordt een klein segmentje van de balk bekeken. In de figuur 3.2 is M(x,t) het moment en stelt V(x,t) de schuifkracht voor.
8
fig. 3.2: Segmentje van de éénzijdig ingeklemde balk. Het verband tussen vervorming en het moment wordt als volgt weergegeven:
EI
d²y = −M dx ²
(3.1)
Veronderstel EI constant, wordt dit gedifferentieerd dan wordt:
EI EI
d³y = −V met V de schuifkracht dx ³
d4y = −W dx 4
met W de lading per eenheid lengte
(3.2) (3.3)
Voor vrije transversale trillingen is de intensiteit van de lading gelijk aan de traagheidskracht. Vergelijking 3.3 wordt dus: d4y d²y EI 4 = − ρA dx dt ²
(3.4)
Dit wordt vereenvoudigd tot: a²
d4y d²y EI + = 0 met a ² = 4 dx dt ² ρA
(3.5)
veronderstel dat y = X(x).T(t)
(3.6)
9
dan wordt dit: (3.7)
....
XT&& + a ² X T = 0 of ....
(3.8)
X T&& =− X a ²T
De linkerkant staat in functie van X terwijl de andere kant in functie van T staat. Iedere kant wordt gelijk gelaten aan de constante ω²/a². Dit geeft volgende vergelijkingen tot gevolg: X − (ω ² / a ²) X = 0
(3.9)
T&& + ω ²T = 0
(3.10)
ω/a= k²
(3.11)
....
De oplossing hiervan is: X(x)= C1cos kx + C2 sin kx + C3 cosh kx + C4 sinh kx
(3.12)
Randvoorwaarden leggen relaties vast met betrekking tot de constanten. Het gaat hier over een aan één zijde ingeklemde balk. Nu is x = 0; X (0) = 0; X& (0) = 0 . M.o.: cosh x =
e x + e−x 2
e x − e−x 2
(3.13)
d (sinh x) = cosh x . dx
(3.14)
sinh x =
en
d (cosh x) = sinh x en dx
Met x=0 in vgl. 3.12 wordt dit C1 + C3 = 0
(3.15)
Wordt vgl. 3.12 gedifferentieerd en x door 0 vervangen, dan wordt dit: C2 + C4 = 0
(3.16)
10
Bij x = L is het moment = 0.Wordt nu vgl. 3.12 twee keer gedifferentieerd, k² verwaarloosd en x vervangen door l dan wordt volgende vergelijking bekomen: 0= -C1cos kl - C2 sin kl + C3 cosh kl + C4 sinh kl
(3.17)
Wordt nu vgl.3.17 drie keer gedifferentieerd, k³ verwaarloosd en x vervangen door l, dan wordt bekomen: 0= C1sin kl - C2 cos kl + C3 sinhkl + C4 cosh kl
(3.18)
Nu is –c1=c3 en c2= -c4 dan worden de 2 laatste vergelijkingen:
Er is een oplossing als de determinant gelijk is aan 0.
Dit zal kloppen voor volgende waarden van kl: 1.875,4.644,7.854,…
Wetende dat:
ω= ak² [rad/s] a² =
(3.19)
EI en kl = 1.875 ρA
(3.20)
Bij het uitrekenen wordt er ook nog omgezet van rad/s naar Hz, dus delen door 2π. 1
1.875² EI 2 F1 = ( ) 2πL ² ρA
(3.21)
Met de gekende parameters kunnen er dus voor de balk verschillende eigenfrequenties uitgerekend worden. Zoals reeds in hoofdstuk 2 vermeld is, kan dit niet toegepast worden op de skilat. Doordat de ski uit verschillende materialen opgebouwd is, moeten de eigenfrequenties experimenteel bekomen worden.
11
4 HOOFDSTUK 4: VOORSTELLEN VAN SYSTEMEN[7] Het proces dat geregeld moet worden, wordt voorgesteld aan de hand van een wiskundig model of transferfunctie. De transferfunctie is het verband tussen ingang en uitgang.
fig. 4.1: Algemene voorstelling van een transferfunctie van een proces. De regeltechniek maakt gebruik van lineaire modellen. Niet-lineaire modellen en processen dienen dus gelineariseerd te worden.
4.1
Voorstelling in het s-domein 4.1.1 Laplace transformatie.
De veranderingen van een proces worden weergegeven door een differentiaalvergelijking, dus in het t-domein. Met de hulp van de Laplace transformatie is het de bedoeling om de moeilijke differentiaalvergelijking te vervangen door een algebraïsche vergelijking. ∞
F ( s ) = ∫ f (t )e − st dt = L{ f (t )}
(4.1)
0
praktisch wordt dit:
d s= dt
t
en
1 = dt s ∫0
(4.2)
Dit geeft als resultaat een transferfunctie van een systeem, en op die manier wordt het dynamisch gedrag van het systeem beschreven. G(s) =
output Y ( s ) = input R( s)
(3.3)
Dit is wel enkel van toepassing op lineaire systemen.
12
4.1.2 Betekenis van nullen en polen De polen zijn de nulpunten van de noemer. Ze bepalen de aard van de tijdsresponsie. De nullen geven het gewicht ervan weer en zijn de nulpunten van de teller. Met de polen kan nagegaan worden of een systeem stabiel is of niet.
Het nulpunten-polen-diagram geeft de nulpunten en de polen van het systeem grafisch weer in het complex vlak. Een nulpunt wordt aangeduid met een ‘o’ en een pool met een ‘x’. Elke pool a komt overeen met een exponentieel tijdsgedrag:
1 → e at s−a •
(4.4)
Indien a positief is, dan wordt de exponentiële functie steeds groter. Het transiënt gedrag sterft niet uit of het systeem is instabiel.
•
Indien a negatief is, dan sterft de reactie vanzelf uit zodat het systeem stabiel is. (Met reactie wordt bedoeld: de uitgang als respons op een willekeurige ingang, bijvoorbeeld ruis). Naarmate a meer negatief is, zal het overgangsverschijnsel vlugger uitsterven.
•
Indien a = 0, dan bevindt het systeem zich op de rand van de stabiliteit. Het overgangsverschijnsel sterft niet uit en wordt ook niet noodzakelijk oneindig groot.
De ligging van de polen en de nulpunten geeft dus informatie over de vorm van de transferfunctie en over het tijdgedrag van het systeem. Figuur 4.2 geeft het nulpuntenpolen-diagram van het eerste orde systeem weer.
fig. 4.2: Nulpunten-polen-diagram van een 1e orde systeem.
13
Onderstaande figuur geeft het verband tussen de ligging van de pool en het transiënt gedrag.
fig. 4.3: Verband tussen de ligging van de pool en het transiënt gedrag. Het proces is: •
Stabiel zolang de polen in het linkerhalfvlak liggen
•
Marginaal stabiel indien de polen zich op de imaginaire as bevinden
•
Instabiel vanaf het moment dat er een pool zich in het rechterhalfvlak bevindt.
4.1.3 Voorstellen in het frequentiedomein. In rusttoestand genereren sinusvormige ingangssignalen bij een lineair systeem sinusvormige responsies. Deze responsies hebben dezelfde frequentie als het ingangssignaal. Ze verschillen echter in fasehoek en amplitude. Dit wordt weergegeven in een bodekarakteristiek. De wijziging van amplitude wordt versterking genoemd. In een bodediagram wordt de as geschaald op 20 log K. De omzetting van een transferfunctie in het Laplacedomein naar het frequentiedomein kan met de vergelijking: s = jω
(4.5)
Een transferfunctie wordt in het frequentiedomein:
G ( jω ) = R(ω ) + jX (ω )
(4.6)
G( jω ) = R²(ω ) + X ²(ω )
(4.7)
De versterking is dan:
14
En de fasehoek:
ϕ = bgtg
X (ω ) R(ω )
(4.8)
fig. 4.4: Algemene voorstelling van een bodekarakteristiek. Een bodediagram kan aan de hand van een rekenprogramma voor elke frequentie berekend worden. Onderstaande figuur geeft de relatie weer tussen de tijdsgrafiek, nulpunten-polen-diagram en het bodediagram.
fig. 4.5: Relatie tussen tijdsgrafiek, nulpunten-polen-diagram en de bodediagram. Ook uit de bodekarakteristiek kan de stabiliteit van het systeem bepaald worden. Er wordt daarvoor gebruik gemaakt van:
15
•
De amplitudemarge AM geeft weer hoeveel de versterking mag toenemen opdat het systeem nog stabiel zou zijn. Hoe groter de amplitudemarge, hoe stabieler het systeem. Het systeem is marginaal stabiel indien de amplitudemarge nul is.
•
De fasemarge FM geeft weer hoeveel de fase mag toenemen opdat het systeem nog stabiel zou zijn. Hoe groter de fasemarge, hoe stabieler het systeem. Het systeem is marginaal stabiel indien de fasemarge nul is.
fig. 4.6: AM en FM aflezen uit een bodekarakteristiek.
4.2
Voorstelling in het Z-domein 4.2.1 Algemeen
De in- of uitgangsgrootheden van een analoog systeem kunnen in principe binnen bepaalde grenzen alle mogelijke waarden aannemen met een oneindige resolutie. Bovendien kunnen deze waarden op elk mogelijk tijdstip veranderen. Men spreekt dan ook van analoge grootheden of signalen. Bij een digitaal systeem kunnen de grootheden slechts op welbepaalde tijdstippen veranderen en heeft de waarde van de grootheid slechts een beperkte resolutie. De overgang van analoge naar digitale signalen noemt men bemonsteren. De tijd tussen twee discrete waarden is de bemonsteringstijd of
16
bemonsteringsperiode, ook wordt dit de sampletijd Ts genoemd. De totale tijd van het begin tot het einde wordt de bemonsteringsduur genoemd.
fig. 4.7: Overgang van analoog naar digitaal signaal.
Hoe kleiner de Ts , hoe beter het analoog signaal benaderd wordt. Er moet wel opgelet worden, bij een te kleine Ts kunnen er problemen optreden met het algoritme.
4.2.2 De Z – transformatie De z-operator vervult bij digitale systemen een analoge rol als de s-operator bij analoge systemen. Om tot de definitie van de Z-getransformeerde te komen, wordt er een ideale bemonsteraar beschouwd. Dan laat het bemonsterde signaal f*(t) zich als volgt schrijven: ∞
f * (t ) = ∑ f (nT )δ (t − nT )
(4.9)
n=0
met:
δ (t ) = dirac impuls
Wordt de laplace-getransformeerde van vgl.4.9 berekend, dan wordt bekomen: ∞
L[ f * (t )] = F ( z ) = ∑ f (nT )e − nsT
(4.10)
n =0
De z-operator wordt als volgt gedefinieerd:
z = e sT
(4.11)
Dit geeft als z-getransformeerde van f*(t): ∞
Z [ f * (t )] = F ( z ) = ∑ f (nT ) z − n
(4.12)
n =0
17
Hieruit kan besloten worden dat de z-getransformeerde van een functie een oneindig lange rij is, maar voor vele standfuncties zijn gesloten uitdrukkingen in z mogelijk.
Stel als functie een éénheidsstap f(t)=1 dan is ∞
F ( z ) = ∑1z − n = 1 + n =0
1 1 + + .... z z²
(4.13)
Deze reeks is convergent voor
z >1 en kan dus in gesloten vorm geschreven worden als
F ( z) =
z z −1
(4.14)
4.2.3 Betekenis van nullen en polen De betekenis ervan is net hetzelfde als bij een analoog systeem. Alleen worden de polen nu uitgezet t.o.v. de eenheidscirkel. De figuur hieronder geeft hierover duidelijkheid.
fig. 4.8: Relatie tussen het Z-domein en het S-domein. Wanneer nu de polen in de cirkel liggen kan er gezegd worden dat het systeem stabiel is. Liggen er polen op de cirkel,dan is het marginaal stabiel en erbuiten onstabiel.
4.2.4 Blokschema van een digitaal regelschema. Een algemeen schema kan eruit zien als de figuur 4.9 . De regelkring opgebouwd voor de actieve demping van de ski is hierop gebaseerd. De computer, met intern de regelaar zal een digitaal signaal uitsturen dat geconverteerd wordt naar een analoog signaal zodat de actuator kan aangestuurd worden. De sensor zal zijn versnelling analoog weergeven,
18
daarom moet dit nog omgezet worden naar een digitaal signaal zodat de PC hier ook mee kan werken.
fig. 4.9: Blokschema van een digitale regelaar.
Digitale regelsystemen worden vooral gebruikt omdat •
er een grote flexibiliteit is t.o.v. analoge regelaars.
•
er weinig problemen zijn met veroudering van componenten en temperatuursinvloed.
•
ze zorgen voor een gemakkelijke interfacing.
•
er gemakkelijk data kunnen gelogd worden.
•
ze minder stoorgevoelig zijn.
19
5 HOOFDSTUK 5: SYSTEEM IDENTIFICATIE, THEORETISCHE BENADERING.[8] 5.1
Algemene theorie.
De system identification toolbox van Matlab/Simulink® laat toe een model op te stellen, gebaseerd op de gemeten input en output data. Hierbij worden nu de parameters aangepast van het model, tot de geschatte output ervan overeenstemt met de gemeten output. Daarvoor kan gebruik gemaakt worden van de graphical user interface, die de meeste functies van de toolbox weergeeft. Zo wordt er vlot toegang verkregen tot alle variabelen die gebruikt worden gedurende het proces. Een andere manier is zelf de commando’s invoeren.
5.1.1 De signalen. Een model zal de relatie tussen 2 gemeten signalen beschrijven, de output signalen worden dan gevormd in functie van de input. Meestal wordt het uitgangssignaal nog beïnvloed door andere niet opgemeten signalen, dit wordt dan ruis of storing genoemd. De onderstaande figuur geeft dit weer.
fig. 5.1: Algemene signaalvoorstelling met u=input,y=output en e=ruis. Deze signalen staan in functie van de tijd, op een tijdstip t worden de waardes gelijk aan u(t), e(t) en y(t). Meestal worden alleen de discrete tijdspunten opgemeten, met als interval de bemonsteringstijd. Het model moet de relatie tussen deze 3 signalen weergeven. Lineaire modellen kunnen als volgt voorgesteld worden in het tijdsdomein: y(t) = Gu(t) + He(t)
(5.1)
20
G duidt de dynamica van het systeem aan, hoe de output gevormd wordt uit de input. Bij lineaire systemen wordt dit de transferfunctie genoemd. H verwijst naar de invloed van de ruis op de uitgang.
5.1.2 Polynoomvoorstelling van een transferfunctie. De functie G en H kunnen ook in termen van z-n weergegeven worden. Modellen berekend via de ARX methode worden zo voorgesteld. G (q ) = q − nk H (q) =
B(q) A(q )
1 A(q )
(5.2) (5.3)
Met A en B als polynomen: A(z)= 1 + a1z-1 + … + anaq-na
(5.4)
B(z)= 1 + b1z-1 + … + bnbq-nb
(5.5)
De termen na en nb duiden op het order van de polynomen waar de term nk het getal is dat het aantal vertragingen van in- tot uitgang weergeeft. Bij het opstellen van een model wordt uitgegaan van de volgende vergelijking: A(z) y(t) = B(z)/F(z) u(t-nk) + C(z)/D(z) e(t)
(5.6)
Dit wordt in de volgende figuur weergegeven :
fig. 5.2: Algemene polynoomvoorstelling A(z),B(z),C(z),D(z) en F(z) zijn hier dus de polynomen. Bij de verschillende modellen die gebruik maken van de polynoomvoorstelling zullen maar enkele van de polynomen aanwezig zijn.
21
Bij de ARX nc,nd en nf gelijk aan 0. Bij de ARMAX nd en nf gelijk aan 0. Bij de ARARX nc en nf gelijk aan 0. Bij ARARMAX nf gelijk aan 0. Bij output-error na,nc en nd gelijk aan 0. Bij Box-jenkins na gelijk aan 0.
5.1.3 Toestandsruimtevoorstelling van een transferfunctie. Een veel voorkomende voorstelling van een lineair systeem is de toestandsruimtevoorstelling: x(t+1) = Ax(t) + Bu(t) + Ke(t)
(5.7)
y(t)= Cx(t) + Du(t) + e(t)
(5.8)
De nx-gedimensioneerde toestandsvector x(t) definieert de relatie tussen u(t) en y(t). De transferfunctie wordt dan: G(z) = C(zInx – A)-1 B + D
(5.9)
G(z) wordt uitgedrukt in functie van de matrices A,B,C en D. H(z) = C(zInx – A)-1 K + Iny
5.2
(5.10)
Modellen voor het identificatieproces.
Het is de bedoeling om het proces dus voor te stellen met een wiskundig model dat zo goed mogelijk het werkelijke proces benadert. Er zijn meerdere modellen , maar hier worden alleen de twee belangrijkste besproken, nl. ARX en N4SID.
5.2.1 ARX Deze methode wordt zeer veel toegepast en blinkt uit door zijn eenvoud. Het nadeel ervan is het feit dat de ruis deel uitmaakt van de systeemdynamica.
22
Er moet wel rekening mee gehouden worden dat deze methode een discreet model weergeeft.
Vertrekkende van de algemene formule wordt dit voor ARX: A(z) y(t) = B(z) u(t-nk) + e(t)
(5.11)
G(z) = z-nk B(z)/A(z)
(5.12)
H(z) = 1/A(z)
(5.13)
Dan wordt :
fig. 5.3: ARX voorstelling.
De meest gebruikte modelstructuur voor een eenvoudige lineaire vergelijking is: y t ( ) a1y + y (t – 1)+…+anay(t – na)=b1u(t – nk)+…+bnbu(t – nk–nb+1)
(5.14)
Deze brengt de huidige output y(t) in relatie tot een oneindig aantal outputs van het verleden y(t-k) en inputs u(t-k). De structuur wordt dus volledig bepaald door na, nb en nk. Na staat voor het aantal polen, nb-1 voor het aantal nullen, terwijl nk de time-delay weergeeft van het systeem. Als er geen dode tijd is wordt nk gelijk aan 1. Om de parameters ai en bi van het arx model te schatten, wordt het volgende commando gegeven: m = arx ( data, [ na nb nk])
Dit geeft de volgende polynomen tot resultaat: A(z)= 1 + a1z-1 + … + anaq-na
(5.15)
B(z)= 1 + b1z-1 + … + bnbq-nb
(5.16)
23
Er zijn nog verschillende meer gespecificeerde modellen, die gedetailleerder zijn dan het basis ARX model. De meest gekende zijn ARMAX, Output-Error en Box-Jenkins.
5.2.2 N4SID Hier wordt gebruik gemaakt van een toestandsruimtevoorstelling om het model te schatten. Dit wordt in deze vorm weergegeven: x(t+1) = Ax(t) + Bu(t) + Ke(t)
(5.17)
y(t)= Cx(t) + Du(t) + e(t)
(5.18)
commando : m = armax ( data, order)
5.3
Graphical User Interface
Met het commando ‘ident’ wordt de graphical user interface(GUI) geopend die de interactie met de toolbox verzorgt. Met dit venster worden de belangrijkste toolbox functies opgeroepen.De gecreëerde variabelen gedurende een sessie kunnen opgeroepen worden via het scherm van de GUI.
24
fig. 5.4: Graphical user interface die de interactie met de toolbox verzorgt. Systeem identificatie is erop gericht om modellen te schatten uit opgemeten data. Daarom is het hoofdvenster van de ident opgesplitst in 2 delen:
Het linkse deel zal het invoeren van data verzorgen Het rechtse deel zorgt voor het verwerken van modellen. De data zullen meestal geïmporteerd worden in de workspace van matlab.
5.4
Datavoorstelling.
In de system identification toolbox worden signalen en data voorgesteld als kolomvectoren.
u (1) u (2) u = ... ... u ( N )
(5.19)
N is het aantal bemonsteringen waarmee de data gemeten zijn, dus bepaald door de bemonsteringtijd Ts. De ingang wordt voorgesteld door u en de uitgang door y.
25
De input-output data wordt voorgesteld door een matrix. De eerste kolom stelt de uitgang voor, de tweede stelt de ingang voor. z=[yu]
(5.20)
Op deze manier kan er makkelijker gewerkt worden met beide signalen.Wanneer de GUI gebruikt wordt voor de identificatie kunnen er stappen gevolgd worden om het resultaat te bekomen. De identificatie kan uitgewerkt worden in 6 stappen[8]:
1. Data importeren. 2. Data bestuderen. 3. Voorbewerken van data. 4. Modellen schatten m.b.v. data. 5. Modellen analyseren. 6. Exporteren van het resulterende model voor verder gebruik.
5.4.1 Data importeren Bovenaan het ident-venster kan er gekozen worden voor data import. Volgende gegevens worden gevraagd: •
Input
•
Output
•
Naam data
•
Starttijd
•
bemonsteringstijd
26
fig. 5.5: Data importeren in de GUI. De input en output data zijn beschikbaar in de workspace van matlab. Opgemeten signalen worden dus eerst in de workspace opgeslagen.Om af te sluiten wordt voor import gekozen.
5.4.2 Data bestuderen Door plot data en data spectra te kiezen kunnen de signalen zowel in het tijds- als in het frequentiedomein bestudeerd worden. Zo kan er bepaald worden of de data voor de identificatie nog moeten worden voorbewerkt.
27
fig. 5.6: Data bestuderen in de GUI.
5.4.3 Voorbewerking van data De gemeten gegevens hebben offset, uitlopers en soms periodes van het missen van waarden. Dit kan tot een incorrect geïdentificeerd systeem leiden. De gemeten gegevens kunnen voorbewerkt worden zodat deze fouten verwijderd of tot het minimum herleid worden. Enkele bewerkingen zijn •
Detrending: Voor het verwijderen van drift en offset.
•
Filteren: om belangrijke frequentiewaaiers in de gegevens te benadrukken en storingen te verwijderen.
•
De sampletijd opnieuw instellen om de schattingstijd en nauwkeurigheid te verhogen.
•
Een bepaald gebied van de metingen selecteren.
5.4.4 Modellen schatten m.b.v. data Hierboven zijn de verschillende modellen reeds besproken. De GUI geeft ons de volgende keuze:
28
fig. 5.7: Gewenst model met zijn orders ingeven. •
Parametrisch model:
Een modelstructuur kenmerkt het verband tussen input en outputgegevens en tussen onbekende ruis en de outputgegevens. Er kan gekozen worden voor: -ARX -ARMAX -Output-Error -Box-Jenkins -State Space •
Spectral model
Geeft meteen een schatting m.b.v. frequentieresponsie. •
Correlation model
Geeft meteen een schatting m.b.v. impulsresponsie.
29
5.4.5 Modellen analyseren De toolbox functies(model view) zorgen ervoor dat de geschatte modeloutput vergeleken kan worden met de output van de data. Zo is zeker dat het geschatte model nauwkeurig het systeem voorstelt. De toolbox verstrekt zes analyses om de geschiktheid van het geïdentificeerde model te bepalen: •
Model output
•
Model residual
•
Transient response
•
Frequency response
•
Zero and poles
•
Noise spectrum
Als er gekozen wordt voor het ARX-model met orders 4 4 1, kan met behulp van model output het gesimuleerde model vergeleken worden met de gemeten output. Door de juiste instellingen aan te brengen bij een bepaald model wordt er een zo groot mogelijke gelijkenis nagestreefd.
fig. 5.8: Gelijkenis controleren tussen model en werkelijke signaal.
30
6
HOOFDSTUK SYSTEEM.[6]
6:
IDENTIFICATIE VAN
2DE
EEN
ORDE
Een 2de orde systeem wordt geïdentificeerd door het opmeten van de in- en uitgang.Zo kunnen de identificatietechnieken gecontroleerd worden alvorens ze toe te passen op de ski. Het 2de orde systeem kan ingesteld worden met behulp van de frequentie F en de demping D. Hierdoor is het systeem gekend en kan het vergeleken worden met het berekende model. De juistheid van de techniek kan dus meteen nagegaan worden.
6.1
Meetopstelling.
1
2 3
4
fig. 6.1: Meetopstelling voor de identificatie van een gekend systeem. 1) Function Generator FG601 2) Scoop Gould Classic 6100 3) 2de orde systeem 4) Voeding
31
6.2
Instelling.
De functiegenerator stuurt een blokgolf uit met een amplitude van 2V en een frequentie van 0.03 Hz. De 2de orde transferfunctie wordt ingesteld met een demping van 0.3 en een frequentie van 0.3. Er wordt een beeld van 50 seconden genomen op de scoop. Met behulp van dSPACE controlldesk® worden de signalen binnengelezen in de pc zodat deze met MATLAB/Simulink® verder verwerkt kunnen worden. Doordat de demping en de frequentie ingesteld is op 0.3, kan nu de transferfunctie berekend worden. Een tweede orde systeem wordt als volgt voorgesteld: G (γω ) =
1
(6.1)
1 + 2 DTγω + T ²ω ²
met T =
1 − D² 2πF
(6.2)
De transferfunctie wordt dus:
TF =
6.3
3.906 s ² + 1.186 s + 3.906
(6.3)
Identificatie.
Zoals reeds in het vorige hoofdstuk besproken, kan de identificatie met verschillende methoden gedaan worden. De identificatie zal m.b.v. de meest gebruikte methoden ARX en n4sid gebeuren. Eerst wordt het bodeplot van het werkelijke systeem vergeleken met het bodeplot van het model. Er moet dus eerst een plot gevonden worden met de in- en output. Dit gebeurt met een in matlab aangemaakt M-file: generatebode.m.(bijlage pagina 15) Om tot een bodekarakteristiek te komen moet er eerst een frequentieanalyse gebeuren. (bijlage pagina 16) Figuur 6.2 toont de bodekarakteristiek van opgemeten data.
32
fig. 6.2: Bodekarakteristiek van het 2de orde systeem. Hierin is een 2de orde systeem te herkennen.Laagfrequent is er een versterking van 0dB en een faseverschuiving van 0°. Hoogfrequent is er een helling van -40dB/decade en een faseverschuiving van -180°. Er is wel veel ruis aanwezig vanaf ± 6 rad/s.
6.3.1 ARX Het 2de orde systeem wordt m.b.v. de ARX-methode geïdentificeerd. De signalen worden i.f.v. de sampletijd in de grafiek uitgezet. Aangezien Ts=0.005 zal 50 seconden gelijk zijn aan 10000 samples. Het programma voor de identificatie is aanwezig in de bijlage.(bijlage pagina 12)
33
fig. 6.3: Voorstellen van ingang en uitgang van het 2de orde systeem. De ingang en de uitgang worden samengevoegd tot een dataobject “ze”. Na het invullen van de orders wordt een “compare” gemaakt tussen het arx-model en de werkelijke uitgang. Als deze gelijkenis voldoende is, wordt er met dit model verder gewerkt. Het geschatte model volgt voor 97.8% het oorspronkelijke signaal.
fig. 6.4: Gelijkenis controleren tussen model en werkelijk signaal.
34
Bij de ARX-methode worden A(q) en B(q) gehaald uit het gekozen model. A(q)y(t) = B(q)u(t) + e(t)
(6.4)
A(q) = 1 - 0.3913 q^-1 - 0.2518 q^-2 - 0.195 q^-3 - 0.1328 q^-4 - 0.134 q^-5 - 0.1023 q^-6 - 0.2868 q^-7 + 0.05544 q^-8 + 0.0635 q^-9 + 0.08631 q^-10 + 0.1528 q^-11 + 0.1373 q^-12
B(q) = 0.0004383 q^-2 - 0.0004443 q^-3 + 0.0006602 q^-4 - 0.0001199 q^-5 + 0.0009429 q^-6 - 0.0009126 q^-7 + 0.0008976 q^-8 - 0.0004646 q^-9 + 0.0003149 q^-10 + 0.0003481 q^-11 - 0.0001119 q^-12 - 0.0002045 q^-13 M.b.v. matlab wordt de transferfunctie berekend:
TF =
− 0.0001659s ² + 0.008304s + 3.558 s ² + 1.359s + 3.573
=
3.558 s ² + 1.359s + 3.573
(6.5)
(6.6)
De gevonden TF is ± gelijk aan de exacte TF, de gebruikte methode heeft tot een juiste oplossing geleid. Hieronder is de bodekarakteristiek van het gevonden ARX-model zichtbaar. Hier is duidelijk een 2de orde systeem zichtbaar, ook is de gelijkenis tussen de opgemeten bode van fig. 6.2 zichtbaar.
35
fig. 6.5: Bodeplot van het 2de orde model.
6.3.2 N4sid Met deze methode worden de resultaten weergeven in een toestandsruimtevoorstelling.
Toestandsmodel: x(t+Ts) = A x(t) + B u(t) + K e(t) y(t) = C x(t) + D u(t) + e(t)
(6.7) (6.8)
De identificatie wordt op dezelfde manier gedaan als bij de ARX-methode. Het programma (bijlage pagina 13) is identiek als dit van de ARX, alleen wordt ARX hier vervangen door n4sid. Na het invullen van de orders wordt een “compare” gemaakt tussen het n4sid-model en de werkelijke uitgang. Als deze gelijkenis voldoende is, kan er met dit model verder gewerkt worden. Het geschatte model volgt voor 98.6% het oorspronkelijke signaal.
36
fig. 6.6: Gelijkenis controleren tussen model en werkelijk signaal. De matrices A, B, C, C, K en x(0) worden in matlab berekend:
M.b.v. matlab wordt de transferfunctie berekend:
TF =
0.0001529s ² + 0.003324s + 3.466 s ² + 1.246s + 3.49
(6.9)
3.466 s ² + 1.246s + 3.49
(6.10)
TF =
Figuur 6.7 toont de bodekarakteristiek van het gevonden n4sid-model. Hier is ook een 2de orde systeem zichtbaar, ook hier is de gelijkenis tussen de opgemeten bode van het werkelijke systeem zichtbaar.
37
fig. 6.7: Bodeplot van 2de orde model.
38
7 HOOFDSTUK 7: IDENTIFICATIE VAN EEN 1E ORDE SYSTEEM. De identificatie van een 1e orde systeem kan op dezelfde manier zoals in het vorige hoofdstuk geïdentificeerd worden. In dit hoofdstuk wordt ook nog een andere manier gebruikt om het beoogde resultaat te verkrijgen. De opstelling is net zoals bij het 2de orde systeem.
7.1
Instelling
De functiegenerator stuurt een blokgolf uit met een amplitude van 2V en een frequentie van 0.03 Hz. De 1e orde transferfunctie wordt ingesteld met een tijdsconstante. Er wordt een beeld van 50 seconden genomen op de scoop. Met behulp van dSPACE worden de signalen binnengelezen zodat deze met matlab verder verwerkt kunnen worden. De tijdsconstante wordt ingesteld op één seconde, zodat de transferfunctie vlot kan berekend worden. Een 1e orde systeem wordt als volgt voorgesteld:
TF =
1 τs + 1
(7.1)
1 s +1
(7.2)
De transferfunctie met τ=1 wordt dus:
TF =
De figuur 7.1 toont zowel het ingangssignaal als het uitgangssignaal.
39
fig. 7.1: Ingang en uitgang, voorstelling van 1e orde systeem.
7.2
Identificatie.
Eerst wordt kort nog eens de identificatie met behulp van de ARX en n4sid methode gedaan. Ten tweede zal een ander manier om te identificeren toegelicht worden. In figuur 7.2 is de bode karakteristiek van de opgemeten data zichtbaar. Er is duidelijk een 1e orde systeem zichtbaar. .Laagfrequent is er een versterking van 0dB en een faseverschuiving van 0°. Hoogfrequent is er een helling van -20dB/decade en een faseverschuiving van -90°.
40
fig. 7.2: Bodeplot van het 1e orde systeem.
7.2.1 ARX Het programma voor de identificatie is net zoals het programma voor de 2de orde, alleen zijn de parameters aangepast. Bij het commando compare wordt er terug een vergelijking gemaakt tussen het werkelijke signaal en het model. Er is een gelijkenis van 97.4% zodat dit zeker meer dan genoeg is om er verder mee te werken. Het programma is te vinden in de bijlagen. (bijlage pagina 10).
41
fig. 7.3: Gelijkenis controleren tussen model en werkelijk signaal. Het programma in matlab rekent de transferfunctie uit:
TF =
0.001908 + 1.004 s + 1.004 =
1 s +1
(7.3) (7.4)
De gevonden TF is ± gelijk aan de ingestelde TF, de gebruikte methode heeft tot een juiste oplossing geleid. Figuur 7.4 toont de bodekarakteristiek van het gevonden ARXmodel. Hier is duidelijk een 1de orde systeem zichtbaar, ook is de gelijkenis tussen de opgemeten bode van fig. 7.2 zichtbaar.
42
Bodekarakteristiek van model:
fig. 7.4: Bodeplot van 1e orde model.
7.2.2 N4SID Het 1e orde systeem wordt ook nog een met de n4sid methode geanalyseerd. Het programma is ook achteraan in de bijlage te vinden.(bijlage pagina 11) Het model komt voor 97.4% overeen met het werkelijke signaal.
43
fig. 7.5: Gelijkenis controleren tussen model en werkelijk signaal. Het toestandsmodel wordt als volgt voorgesteld: dx/dt = A x(t) + B u(t) + K e(t)
(7.5)
y(t) = C x(t) + D u(t) + e(t)
(7.6)
De matrices en de transferfunctie worden terug matlab berekend zodat het volgende resultaat bekomen wordt:
TF = =
1.003 s + 1.002
(7.7)
1 s +1
(7.8)
Wat ook zeer goed het werkelijke systeem benadert.
44
7.2.3 Identificatie d.m.v. GUI Ook met behulp van de GUI kan de transferfunctie gevonden worden. Oproepen gebeurt via het commando “ident”. Via data import wordt volgend venster verkregen:
fig. 7.6: Data importeren via deGUI. De in- en uitgang worden ingevuld alsook de sampletijd. Als dit geïmporteerd is, wordt er teruggegaan naar het ident-venster. Bij estimate wordt vervolgens parametric models gekozen. Bij structure wordt de ARX gekozen en de gewenste orders ingegeven.
fig. 7.7: Parametric models,keuze van het model.
45
Wanneer nu de model output in het ident venster gekozen wordt, zal de gelijkenis tussen het model en de uitgang van het signaal weergegeven worden.
fig. 7.8: Gelijkenis controleren tussen model en werkelijk signaal. Als er dubbel wordt geklikt op ARX441 in het ident venster, verschijnt info van dit model.
fig. 7.9: Info opvragen van gewenste model.
46
Is de gelijkenis groot, dan worden de nulpunten en de polen opgevraagd. Nu wordt “zeros en poles” aangevinkt.
fig. 7.10: Eénheidscirkel van het gekozen model. Aangezien de dominante polen het transiënt gedrag bepalen, zal er nu enkel rekening worden gehouden met de dominante polen en nullen. Er wordt slecht 1 dominante pool waargenomen, nl. deze op de reële as waarvoor geldt z =0.995. Dit is nu een pool van het digitaal systeem, als nu de pool van het analoog systeem gezocht wordt moet er nog een omzetting gebeuren van het z-vlak naar het s-vlak. Met: z = esT met s = a + bj
(7.8)
De dominante pool kan geschreven worden als z = 0.995,vgl. 7.8 wordt: 0.995 = esT ln 0.995 = sT
met T=0.005
(7.9) (7.10)
Dit wordt: s= -1.0025 ≅ -1
(7.11)
Dus wordt dan de transferfunctie:
TF =
1 s +1
(7.12)
47
8 HOOFDSTUK 8: IDENTIFICATIE VAN DE SKILAT 8.1
Systeembeschrijving.
fig. 8.1: Proefopstelling met actuator en sensor. Dit is de meetopstelling, waarbij de ski ter hoogte van de voet ingeklemd is. De actuator wordt aangestuurd door een sinus(ingang) terwijl de versnellingsopnemer de versnelling zal opmeten(uitgang). Er is reeds aangetoond dat de ski vooral zal trillen met zijn eigenfrequentie van ± 12.5 Hz. Daarom wordt de actuator aangestuurd door een sweepsignaal. Dit is een signaal dat start bij een frequentie van 0 Hz en zo stijgt tot een aan te geven frequentie. Er wordt gekozen voor een eindfrequentie van 25Hz. Dit wordt gemeten over een periode van 10 seconden met een sampletijd van 0.001 seconde. Er waren ook problemen met het versnellingssignaal. Dit signaal was namelijk zeer klein en werd via de dSPACE-kaart binnengelezen. Bij dSPACE controlldesk® wordt het analoog signaal omgezet naar bits, maar omdat het signaal te klein was zonder versterking, was het verkregen signaal zeer hoekig en zeker en vast niet mooi sinusoidaal. dSPACE controlldesk® werkt met een signaal tussen 0 en 1, daarvoor heeft het een bepaald aantal bits ter beschikking. Doordat het signaal zonder versterking zeer klein is, wordt er maar een klein aantal bits gebruikt om het signaal weer te geven. Dit probleem wordt opgelost door een versterker tussen het versnellingssignaal en de dSPACE-kaart te plaatsen. Zo wordt het aantal gebruikte bits dus meteen verhoogd. Dat heeft een mooier signaal tot resultaat, waardoor er makkelijker mee verder kan gewerkt worden. Ook wordt er gebruik gemaakt van de ruisfunctie die toelaat om de ruis wat te filteren. Figuur 8.2 toont het
48
aangelegde signaal en het signaal van de versnellingsopnemer. Het ingangssignaal is dus de sweep zoals zichtbaar is op de figuur. Bij het uitgangssignaal is het duidelijk dat de uitwijking veel groter is bij de eigenfrequentie van de skilat. Het valt op dat er quasi geen respons is voor lage frequenties (<7 à 8Hz). Dit komt omdat de versterking nog altijd niet sterk genoeg is voor de zeer kleine uitwijkingen bij de lage frequenties. De apparatuur om nog meer te versterken was niet beschikbaar, zodat er verder moest gewerkt worden met deze resultaten.
fig. 8.2: In- en uitgang van aangelegd sweep-signaal.
8.2
Identificatie
Het was de bedoeling om ons systeem te identificeren m.b.v. een soortgelijk programma dat gebruikt wordt voor de identificatie van het 1e en 2de orde systeem. Deze methode heeft echter nooit tot een valabele oplossing geleid. Een verklaring hiervoor is misschien dat er quasi geen respons is voor lage frequenties (<7 à 8Hz). Gelijk welk algoritme zal het met deze data moeilijk hebben voor lage frequenties. Verder mag je rekentools niet zomaar klakkeloos gebruiken, althans wat er uit komt, mag niet meteen geloofd worden. De ultieme handleiding voor identificatie bestaat niet. Er zijn natuurlijk wel regels die je
49
meestal tot een aanvaardbaar resultaat brengen, maar zoals dat met elke tool is, kan je rare dingen uitkomen en dat kan moeilijk in een formule of algemene handleiding gestopt worden. Daarom is er om tot een oplossing te komen, gebruik gemaakt van de GUI tools.
8.2.1 Identificatie m.b.v. de GUI. Na het openen van de GUI worden de signalen geïmporteerd. De figuur 8.3 toont de opgevraagde time plot waar terug de in- en uitgangssignalen zichtbaar zijn. De sweep bouwt zijn frequentie op tot 25 Hz. Rond zijn eigenfrequentie van 12.5 Hz is de versnelling duidelijk veel groter.
fig. 8.3: In- en uitgang van aangelegd sweep-signaal in de GUI. Er wordt nu een model geschat via estimate, en er wordt dan gekozen voor parametric models. De n4sid methode wordt gebruikt en de orde wordt gelijk gesteld aan 6.
50
fig. 8.4: Parametric models,keuze van gewenst model. Dit wordt nu via estimate geschat. De uitgang van het model volgt voor 80 % de uitgang van het werkelijk systeem. Moest de respons voor het laag frequente signaal wat groter zijn, zou de gelijkenis hoger liggen dan 80 procent. Als de signalen rond de eerste eigenfrequentie vergeleken worden, zijn de 2 signalen bijna gelijk aan elkaar. Daarom wordt er met dit model verder gewerkt.
fig. 8.5: Gelijkenis controleren tussen model en werkelijk signaal.
51
Nu wordt het nulpunt-polen-diagram van dit model bekeken, figuur 8.6 geeft dit weer. Er wordt alleen maar rekening gehouden met de dominante polen en nullen, aangezien deze het transiënt gedrag bepalen.
fig. 8.6: Voorstelling van de éénheidscirkel van het model.
52
fig. 8.7: Eénheidscirkel model met dominante polen en nullen (ingezoomd). Uit het nulpunten-polen-diagram worden de dominante polen en nullen weerhouden. Dominante polen: 0.995 ± 0.08i Dominante nullen: 1.015 en 0.97
Om nu de transferfunctie te vinden moet nu nog de versterkingfactor K gevonden worden. Deze kan experimenteel bekomen worden door gebruik te maken van Matlab/Simulink®. Een nauwkeuriger manier om deze te vinden is door ze uit de bodekarakteristiek te halen. Met het programma bodekarakteristiek.m wordt figuur 8.8 verkregen.
53
fig. 8.8: Bodeplot van de skilat bekomen met sweep-signalen. Uit de bodekarakteristiek kan nogmaals de eigenfrequentie afgeleid worden. Er is een amplitudepiek bij 78.5 rad/s wat overeenstemt met de reeds gevonden 12.5 Hz. De karakteristiek begint ongeveer bij een amplitude van -14.186 dB en aangezien de polen en nullen in het bodediagram starten bij een amplitude nul, is het dus de versterkingsfactor K die verantwoordelijk is voor deze verschuiving. Dan is 20 log K = −14.186
(8.1)
K = 0.1953
(8.2)
en wordt:
De transferfunctie in de zpk vorm:
0.1953( z − 0.9635)( z − 1.015) ( z ² − 1.988 z + 0.9945)
(8.3)
of
TF =
0.1953z ² − 0.3864 + 0.191 ( z ² − 1.988 z + 0.9945)
(8.4)
54
9 HOOFDSTUK 9: ANALYSE VAN EEN TOESTANDSRUIMTE.[9] Een complex systeem kan vele inputs en outputs hebben waarbij de relatie ertussen zeer gecompliceerd kan zijn. Daarom is het van belang om sommige systemen te vereenvoudigen in hun wiskundige uitdrukking. Voor de verwerking met de computer is het nuttig te beschikken over een gestandaardiseerde oplossingsmethode, onafhankelijk van de complexiteit van het systeem. De analysemethode m.b.v. de toestandsruimte is hiervoor de ideale methode. De toestandsruimtevoorstelling is moderner dan de transferfunctie. Lineaire en niet-lineaire systemen, met 1 of meerdere in- en uitgangen (MIMO-systemen) kunnen voorgesteld worden op een uniforme manier. In de toestandsruimte wordt een systeem beschreven door een set van veranderlijken die eigen zijn aan de inwendige toestand van het systeem op een gegeven ogenblik en door de vergelijkingen tussen deze veranderlijken in functie van de tijd. Het aantal veranderlijken dat nodig is om de toestand waarin het systeem zich bevindt eenduidig te bepalen, is gelijk aan de orde van het systeem. De vergelijkingen tussen de verschillende toestandsveranderlijken zijn steeds differentiaalvergelijkingen of differentievergelijkingen van de eerste orde. Het aantal vergelijkingen is weerom gelijk aan de orde van het systeem. In dit hoofdstuk wordt eerst een algemene voorstelling van de toestandsruimte gegeven. Daarna wordt de toestandsterugkoppeling en het maken van een toestandschatter besproken.
9.1
Algemene voorstelling
Een gewoon differentiële vergelijking, die het gedrag van een lineair systeem weergeeft, wordt meestal als volgt weergegeven: dx(t)/dt= Ax(t) + Bu(t)
(9.1)
y = Cx(t) + Du(t)
(9.2)
De eerste vergelijking wordt de toestandsvergelijking genoemd, de tweede vergelijking is de uitgangsvergelijking welke de relatie tussen de uitgangen, de toestandsvariabelen en de ingangen weergeeft. Daarbij is
x de toestandsvector u de ingangsvector y de uitgangsvector
55
De matrices A,B,C en D variëren hier met de tijd.
Een lineaire differentiële vergelijking met constante coëfficiënten dat het gedrag van een lineair systeem dat niet varieert in de tijd(LTI) weergeeft, ziet er dan als volgt uit: dx/dt= Ax+ Bu
(9.3)
y = Cx + Du
(9.4)
Hier zijn A,B,C en D dus constante matrices. De figuur hieronder toont een blokschema dat de toestandsvergelijkingen voorstelt.
fig. 9.1: Blokschema toestandsvergelijking.
9.2
Wiskundige modellen van lineaire systemen.[1]
Matlab beschikt over een commando om een vlotte omschakeling tussen een transferfunctie en een toestandsruimte te maken en omgekeerd.
van TF naar SS. Het commando [A,B,C,D] = tf2ss(num,den) zet de TF van de vorm Y(s)/U(s)= num/den = C(sI – A)-1B + D
(9.5)
om naar de toestandsruimte: dx/dt= Ax+ Bu
(9.6)
y = Cx + Du
(9.7)
56
van SS naar TF Het commando [num,den] = tf2ss (A,B,C,D)
geeft de TF Y(s)/U(s) weer, indien het systeem slechts 1 input en output heeft. Zijn er meerdere in- en uitgangen dan ziet het commando er zo uit. [num,den] = tf2ss (A,B,C,D,iu)
Dit zet de toestandsruimte dx/dt= Ax+ Bu
(9.8)
y = Cx + Du
(9.9)
om naar de transferfunctie: Y(s)/Ui(s)= num/den = C(sI – A)-1B + D
(9.10)
Met de i wordt er aangeduid welke ingang er gebruikt wordt voor responsie van de uitgang.
9.3
Controleerbaarheid[2]
Elk systeem is controleerbaar op een tijdstip t0 als het mogelijk is om het systeem van zijn begintoestand x(t0) naar een andere toestand te brengen binnen een eindig tijdsinterval. Vertrekkende van de toestandsvergelijking van een nde orde systeem dx/dt= Ax + Bu
(9.11)
Eén van de verschillende methoden om de controleerbaarheid te testen, is het gebruik van een controleerbaarheidsmatrix S = [ B AB A²B A³B … An-1B]
(9.12)
waarvan de rank = n.
In matlab wordt de rank met het volgende commando opgevraagd:
rank(ctrb(A,B))
57
Is dit gelijk aan n: – via de input zijn alle n toestanden stuurbaar. – volledige toestandsterugkoppeling zal zinvol zijn.
9.4
Observeerbaarheid[2]
Een systeem is observeerbaar als er een toestand x(t) kan bepaald worden met de output binnen een eindig tijdsinterval.Vertrekkende van de toestandsvergelijking van een nde orde systeem wordt er nu ook nog de uitgangsvergelijking bijgenomen: dx/dt= Ax+ Bu
(9.13)
y = Cx + Du
(9.14)
Eén van de verschillende methoden om de controleerbaarheid te testen, is het gebruik van een observeerbaarheidsmatrix die een rank = n moet hebben bij D = 0
C CA V = CA² M CA n −1
(9.15)
In matlab wordt de rank met het volgende commando opgevraagd: rank(obsv(A,C)) Is dit gelijk aan n: – alle n toestanden werken door in de output y. – vanuit y is reconstructie van de toestanden, via een schatter, mogelijk.
9.5
Stabiliteit.[7]
Hiervoor wordt terug vertrokken van de toestandsvergelijking: dx/dt= Ax + Bu
(9.16)
het ingangssignaal kan weggelaten worden: dx/dt= Ax
(9.17)
58
de oplossing is van de exponentiële vorm: x = keλt
(9.18)
λkeλt =Akeλt
(9.19)
λx = Ax
(9.20)
(λI-A)x = 0 x = (λI-A)-1 = 0 Dit is de karakteristieke vergelijking.
(9.21)
Dan wordt:
Deze vgl. heeft enkel oplossing als (λI-A) = 0. De oplossing levert de eigenwaarde of de karakteristieke polen op.( λ = de karakteristieke wortels of eigenwaarden)
9.6
Poolplaatsing door toestandsterugkoppeling.
Als het systeem volledig controleerbaar is met de toestandsruimte, dan mogen de polen van het gesloten lus systeem, op een zelf te kiezen plaats, gekozen worden. Dit gebeurt dan via de versterkingsmatrix van de toestandsterugkoppeling. Er wordt terug vertrokken van de volgende toestandsvergelijking: dx/dt= Ax + Bu
(9.22)
u = -Kx + r
(9.23)
Uit de figuur 9.2 blijkt dat:
waarbij r de referentie-input voorstelt, K is de 1x n terugkoppelingsmatrix. K= [ k1 k2 k3 … kn ]
(9.24)
Wordt vgl. 9.23 nu in vgl. 9.22 gestoken dan is het resultaat: dx/dt= Ax + B(-Kx + r)
(9.25)
dx/dt= (A – BK)x + Br
(9.26)
Veronderstel r = 0 dan wordt:
dx/dt= (A – BK)x
(9.27)
eigenwaarden: det(λI - ( A – BK))
59
De matrix ( A –BK) is de systeemmatrix van het gesloten lus systeem. De eigenwaarden ervan worden de regulator poles genoemd en bepalen de dynamica van het systeem. Liggen deze polen in het linker-halfvlak, dan zal x(t) naderen naar 0 terwijl t naar het oneindige evolueert.
fig. 9.2: Toestandsterugkoppeling. Bij het verplaatsen van de polen geldt het volgende: •
Bij het verplaatsen van de polen naar links, wordt het systeem sneller.
•
Bij het verplaatsen van de polen naar een lagere imaginaire waarde, wordt de demping groter.
60
fig. 9.3: Vlugge demping met behulp van het nulpunten-polen-diagram. Is er een snel gedempt systeem gewenst, dan moeten polen naar elkaar toe geschoven worden en naar links geplaatst worden. De rechte lijnen stellen de dempingslijnen voor, de ellipsen geven de pulsatielijnen weer. Het bekomen van de K matrix kan op 2 manieren gebeuren, manueel uitrekenen ofwel matlab het werk laten doen. Dit wordt aangetoond met behulp van een voorbeeld.
2 − 2 A= B= 1 0
1 0
C = [0
1] en D = 0
Er mag aangenomen worden dat dit systeem observeerbaar, controleerbaar en stabiel is. De polen moeten bijvoorbeeld op deze locatie gekregen worden: -5± i Dan is de karakteristieke vgl.: (s + 5+ i)(s + 5- i) = s² + 10 +26 = 0 Wordt nu de KV bepaald van het systeem:
61
det(λI - ( A – BK)) uitrekenen:
2 − 2 1 A – BK = - [k1 1 0 0
2 − 2 k1 k2 ] = - 1 0 0 2 − k1 = 1
k2 0
− 2 − k2 0
dit wordt nu:
λ 0 2 − k1 0 λ - 1
− 2 − k 2 λ − 2 + k1 = 0 − 1
2 + k2 λ
λ+(k1 –2) λ+k2+2 = 0 de 2 KV moeten aan elkaar gelijk zijn λ+(k1 –2) λ+k2+2 = s² + 10 +26 = 0
k1 –2 = 10 k2+ 2 = 26 dus
k1 = 12 en k2 = 24
Dit kan nu ook allemaal via matlab gebeuren, gewoon de matrices invoegen en dan het commando K=place(A,B,[gewenste poollocaties])
K = place(A,B,[-5-i;-5+i])
K=
12 en 24
62
9.7
Toestandschatter.
Indien de toestanden van het systeem niet gemeten kunnen worden, kan er ook geen toestandsterugkoppeling gebeuren. In dit geval moeten de toestanden geschat worden. Dit gebeurt met een schatter die het werkelijke systeem simuleert. De schatter moet zo goed mogelijk overeenstemmen met het werkelijke systeem. De matrices A, B, C en D van de schatter moeten dus overeenstemmen met deze van het systeem. De ingang van de schatter is gelijk aan de ingang van het systeem. Indien de schatter goed werkt moet ook de uitgang van de schatter gelijk zijn aan deze van het systeem. Vaak zullen deze uitgangen echter nog verschillen. Het foutsignaal op de uitgang kan dan gebruikt worden om de schatter te verbeteren en aan te passen. Dit gebeurt via de terugkoppelmatrix L. Het volledige simulatieschema van systeem en schatter is gegeven in de figuur 9.4. Hierbij wordt de matrix D gelijk gesteld aan 0.
fig. 9.4: Algemene voorstelling van toestandsregelaar met toestandsschatting. De schattingsfout wordt gegeven door: ~ x (t ) = x(t ) − xˆ (t )
(9.28)
Wordt nu terug de basis toestandsvergelijking beschouwd: x ( t + 1 ) = Ax ( t ) + Bu ( t )
(9.29)
63
Uit de figuur kan duidelijk nagegaan worden dat:
xˆ (t + 1) = Axˆ (t ) + Bu (t ) + L( y (t ) − yˆ (t )) en
yˆ (t ) = Cxˆ (k )
Het verschil tussen formule 9.29 en formule 9.30 geeft volgend resultaat: ~ x (t − 1) = A( x(t ) − xˆ (t )) − L( y (t ) − yˆ (t ))
~ x (t − 1) = A( ~ x ) − L(Cx (t ) − Cxˆ (t ))
(9.30)
(9.31) (9.32)
x (t − 1) = A~ x (t ) − LC~ x (t ) ~
(9.33)
~ x (t − 1) = ( A − LC ) ~ x (t )
(9.34)
Dit wordt dan:
De terugkoppeling via L moet de fout op de uitgang van de schatter nul maken. Dan kan ervan uitgegaan worden dat ook de toestanden van de schatter de werkelijke toestanden van het systeem weergeven. Deze geschatte toestanden worden teruggekoppeld naar het systeem via K. De gesloten-systeemmatrix blijft dan (A - BK). De evolutie van de fout op de uitgangen wordt beschreven door de matrix A-LC.Op dezelfde manier waarmee de terugkoppelingsmatrix bekomen wordt, kan L bepaald worden door poolplaatsing. De schattingsfout moet nul worden voor t naar ∞, dan moet: det(λI - ( A – LC)) = 0
(9.35)
Dit kan opnieuw makkelijk verkregen worden m.b.v. matlab:
L=place(A,C,Pobsv)
Bij de keuze van de polen van de observer is de vuistregel dat deze 5 à 10 maal meer naar links moeten liggen. Zo wordt het schattingsproces vlugger dan het regelproces, zodat de dynamica van het gesloten lus systeem nauwelijks beïnvloed wordt door de dynamica van de schatter.
64
10 HOOFSTUK 10: ONTWERP TOESTANDSRUIMTEREGELAAR.[7] De toestandsregelaar is een zeer krachtige en flexibele regelaar. Deze moet gezien worden als de tegenhanger van de klassieke PID-regelaar. In tegenstelling tot deze laatste regelaar waarmee enkel SISO systemen geregeld kunnen worden, is de toestandsregelaar ideaal geschikt om ook systemen met meerdere in- en uitgangen te regelen (MIMO systemen).De toestandsregelaar kan als een discreet werkende regelaar of als een continue regelaar ontworpen worden. Vaak zal dit soort regelaars echter m.b.v. een computersysteem geïmplementeerd worden zodat het noodzakelijkerwijs een digitaal werkende regelaar wordt. In feite is de toestandsregelaar niets anders dan een terugkoppeling van de toestanden van het gegeven systeem met de juiste gewichtsfactoren. Hiervoor moet het systeem eerst beschreven worden in de toestandsruimte. Dit resulteert in het toestandsruimtemodel dat zoals de TF een beknopte wiskundige weergave is van het systeem. De volgende paragrafen behandelen het ontwerp van een toestandsruimteregelaar. Er wordt getoond hoe een model vanuit de TF kan opgesteld worden. Vervolgens wordt de invloed van een toestandsterugkoppeling besproken. Door middel van poolplaatsing zal een toestandsterugkoppelmatrix ontworpen worden en zo het volledige systeem naar de hand worden gezet.Tot slot wordt een model gegeven van een toestandsschatter die de toestanden van het systeem zal berekenen indien deze niet gemeten kunnen worden.
In hoofdstuk 8 is de TF voor ons systeem bepaald.
TF =
0.1953z ² − 0.3864 z + 0.191 z ² − 1.988 z + 0.9945
(10.1)
Dit wordt omgezet naar een analoog systeem.
TF =
0.1953s ² + 3.629s − 100.3 s ² − 5.515s + 6521
(10.2)
Dit wordt nu omgezet naar een toestandsruimtemodel. Daarvoor moeten zowel noemer als teller als matrix ingeven worden. Commando’s: b=[0.1953 3.629 -100.3] a=[1 5.515 6521]
65
Met het volgende commando worden de matrices A, B, C en D bekomen. [A,B,C,D]=tf2ss(b,a)
A= 1.0e+003 * -0.0055 -6.5210 0.0010
0
B= 1 0
C= 1.0e+003 * 0.0026 -1.3739
D= 0.1953
Hiermee is het toestandsmodel gekend. Er kan makkelijk geverifieerd worden of het model en de TF(zowel discreet als analoog) gelijk zijn.Er wordt een puls aangelegd en op de scoop kan nagegaan worden of ze alle drie op dezelfde manier reageren.
fig. 10.1: Controle toestandsmodel.
66
fig. 10.2: Scoopbeeld van controle toestandsmodel. Op het scoopbeeld is te zien dat de 3 signalen op elkaar liggen. Er is dus zekerheid dat het berekende toestandsmodel klopt. De waarden van het toestandsmodel worden in Matlab/Simulink® geladen via het M-file toestandsmodel.m.(Bijlage pagina 14)
10.1
Toestandsterugkoppeling.
Vooraleer de terugkoppeling gemaakt wordt, moeten eerst de polen van het systeem van de ski berekend worden. De polen van het gesloten systeem zijn de eigenwaarden van A. Eig(A) = -2.7500 +80.7059i -2.7500 -80.7059i
67
fig. 10.3: Nulpunten-polen-diagram van de skilat. Er wordt een puls van 1 seconde aan het toestandsmodel aangelegd. Het uitgangssignaal zal dus het signaal zijn dat de versnellingsopnemer in het systeem zou opmeten. Zo wordt volgende schema opgebouwd:
fig. 10.4: Schema toestandsterugkoppeling.
68
Merk op dat bovenaan het schema gewoon het systeem wordt geplaatst, zo wordt het makkelijker om te vergelijken tussen het signaal waar poolplaatsing is op toegepast en het gewone signaal. De polen van het systeem zijn gekend. Het nieuwe systeemgedrag wordt bepaald door de nieuwe polen. Eerst wordt de imaginaire waarde niet aangepast, alleen het reële deel wordt veranderd. De polen worden naar links opschoven zodat het systeem veel sneller wordt. De gekozen polen zijn: -20 +80.7i -20 -80.7i
Het scoopbeeld ziet er als volgt uit(geel signaal is niet gedempt, blauwe wel):
fig. 10.5: Scoop 1: Voorstelling van een gedempt systeem. Het systeem is sneller geworden, de pulsaties van de 2 signalen zijn gelijk, wat volledig overeenstemt met hetgeen er gevraagd is. Dit komt omdat we de imaginaire factor niet
69
aangepast hebben. Vervolgens wordt deze factor ook aangepast, ook de reële factor wordt nog iets negatiever gemaakt. De gewenste polen zijn nu: -30 +30i -30 -30i Het scoopbeeld ziet er nu als volgt uit:
fig. 10.6: Scoop 1: Voorstelling een snel gedempt systeem. Er is nu meer demping en de signalen hebben ook niet meer dezelfde pulsatie. Door de polen nog wat op te schuiven naar links is het systeem ook nog iets sneller. Door het veranderen van deze factoren kan er dus een zeer snel systeem ontstaan. Als de polen ver naar links geplaatst worden met een lage imaginaire term, kan er dus een zeer snel gedempt systeem verkregen worden. Een systeem met volgende polen is hiervan een voorbeeld: -100 +i -100 –i
70
Het scoopbeeld:
fig. 10.7: Scoop 1: Voorstelling van een zeer snel en gedempt systeem. Het is dus duidelijk dat met behulp van een toestandsterugkoppeling het systeemgedrag veranderd kan worden. De keuze van de polen zal dus bepalend zijn voor demping van de ski. Op deze manier zou dus de mate van demping door de gebruiker kunnen ingesteld worden.
10.2
Toestandsschatter
Nu wordt de schatter opgebouwd die de toestanden van het systeem zal schatten zodat de toestandsterugkoppeling gebruiken kan worden. Het model voor een schatter werd reeds besproken, nu moet alleen nog rekening gehouden worden met de matrix D die nu niet meer gelijk aan 0 is. De polen van de observer worden gekozen zodat hij ±tien keer sneller is dan het proces.
71
Polen: -30+i -30-i
Het schema is hieronder afgebeeld:
fig. 10.8: Toestandschatter. De scoopbeelden vormen een controle, hiermee kan de werking van de observer nagegaan worden. Scoop 1 toon ons het geschatte signaal. Dit kan vergeleken worden met het echte signaal dat hierboven reeds enkele keren is afgebeeld. Maar om zeer mooi te zien of de schatter werkt, moet scoop 2 bekeken worden. Deze beeldt de schattingsfout af. Op het 2de scoopbeeld kan de schattingsfout afgelezen worden. Deze was ogenblikkelijk gelijk aan nul. Er mag dus gesteld worden dat de schatter goed werkt.
72
fig. 10.9: Scoop1: Uitgangssignaal van de toestandsschatter.
73
10.3
Toestandsruimteregelaar[5]
De regelaar in zijn geheel ziet er dan als volgt uit:
fig. 10.10: Totale toestandsruimteregelaar. De polen van de terugkoppeling zijn nog steeds: -100 +i -100 –i
De polen van de observer zijn nog steeds: -30+i -30-i
Scoop 1 geeft het verschil weer tussen het niet-gedempt signaal en het gedempt signaal.
74
Dit wordt hieronder afgebeeld:
fig. 10.11: Scoop1: Demping van signaal.
Het signaal wordt zeer vlug gedempt. De regelaar doet uitstekend zijn werk. In ons schema hierboven stelt het blokje state-space1 ons systeem voor. In de opstelling van de ski wordt dit vervangen door de ski zelf. Het signaal voor het blokje en dus hetgeen gemeten wordt met scoop 2 is hetgeen gestuurd zal worden naar de actuator. Het signaal na het blokje zal vervangen worden door het signaal van het versnellingsopnemertje. Zoals op scoop 1, is zichtbaar dat het signaal dat de versnellingsopnemer uitstuurt naar nul zal gedwongen worden. Dit zal gebeuren door het signaal op scoop 2 dat dan de actuator zal aansturen. In plaats van een puls aan te leggen, zal bij de ski een constante 0 aangelegd worden, want er wordt geen uitwijking gewenst.
75
fig. 10.12: Scoop2: Stuursignaal van de regelaar.
76
11 HOOFDSTUK 11: ACTIEVE DEMPING VAN DE SKI 11.1
Implementatie van de regelaar.
Voor de demping van de ski moet nu de regelaar toegepast worden op de skilat. Het fysisch systeem kan nu ook in de regelaar gebracht worden. De waardes voor de poolplaatsing blijven dezelfde als bij het laatste voorbeeld uit het vorige hoofdstuk.Er wordt als ingangssignaal 0 aangelegd, want er wordt een uitwijking van 0 geëist. 0 is dus de gewenste waarde die aan de regelaar aangelegd wordt. De versnelling wordt ook aangebracht aan de regelaar. Deze zal vergeleken worden met de uitgang van de toestandsschatter. Aangezien deze veel sneller is door de poolplaatsing, zal er dus een fout zijn die de regelaar zo vlug mogelijk moet tegenwerken zodat de fout gelijk wordt aan 0. Daarom zal de actuator met de aanstuurspanning worden aangestuurd om op gepaste wijze tegen te trillen. Dit is de ingang van ons systeem. Figuur 11.1 toont het totale regelschema voor de actieve demping van de skilat.
fig. 11.1: Toestandsruimteregelaar van de skilat.
77
Dit schema kan voorgesteld worden door figuur 11.2. Op deze figuur is de algemene voorstelling van de regelaar in figuur 11.1 zichtbaar.
fig. 11.2: Algemene voorstelling van de regelaar
11.2
Controle werking.
De werking wordt gecontroleerd door een trilling te genereren. In figuur 11.3 worden de trillingen, zowel gedempt als niet-gedempt, op hetzelfde ogenblik opgemeten. Deze figuur toont zeer mooi het verschil tussen de twee trillingen. Er mag gezegd worden dat de regelaar uitstekend werkt.
fig. 11.3: Vergelijking gedempte en niet-gedempte trilling.
78
De goede werking van de regelaar is dus aangetoond. Tot slot moet nog even opgemerkt worden dat actieve demping van een ski natuurlijk niet met behulp van een trillingsopwekker kan gebeuren. Bij ski’s zal men piezo-elektrisch materiaal inbrengen in de ski om alle trillingen weg te werken en zo het contact met de sneeuw te behouden. Hierbij zal het piëzo-elektrisch materiaal elke trilling omzetten in een spanning via een elektronisch circuit. Dit circuit stuurt dan een tegenspanning naar het piëzo-elektrisch materiaal wat een tegengestelde kracht opwekt en zo de trillingen uitdempt.
79
BESLUIT. Het doel van dit eindwerk bestaat erin om ongewenste trillingen zo vlug mogelijk uit te dempen. Hiervoor is een experimentele proefopstelling gebouwd m.b.v. een skilat. De opstelling bestaat uit een skilat met daarop een versnellingsopnemer gemonteerd, die de versnelling bij een trilling weergeeft. Er wordt ook een actuator aan de skilat bevestigd. In eerste instantie bevat dit eindwerk een theoretische studie van de verschillende identificatietechnieken die MATLAB® aanbiedt. Vooraleer de trillingen gedempt kunnen worden, is het namelijk noodzakelijk om een gepast model op te bouwen van ons systeem. Het is belangrijk om te weten hoe de skilat reageert op verschillende trillingen. Dit gebeurt door het opmeten van de responsie wanneer er geëxciteerd wordt met een multisinus. Het excitatiegebied is 25 Hz, zodat alleen de eerste eigenmode van 12.5 Hz hierbinnen zal vallen. Met behulp van dit nauwkeurig wiskundig model, is het mogelijk een toestandsruimteregelaar te ontwerpen. De terugkoppeling wordt ontwerpen met de poolplaatsingstechniek. Hierbij worden de polen zo gekozen dat er een demping van de eerste eigenmode optreedt. De toestanden zijn niet fysisch meetbaar, daarom wordt een toestandsschatter gebruikt. Deze schat de toestanden van het systeem op basis van het gevonden wiskundig model. Dit alles zit omvat in de regelaar. De ongewenste trillingen van de skilat worden op deze manier gedempt. Zowel de amplitude als de duur van een trilling wordt op deze manier serieus verminderd. Dit eindwerk heeft het dempen van een skilat uitvoerig besproken. Bij andere materialen kan dit evenzeer toegepast worden. Het is wel uiterst belangrijk om het te dempen materiaal eerst aan een serieuze analyse te onderwerpen zodat het wiskundig model zo goed mogelijk het werkelijke systeem benadert.
80
Literatuurlijst. [1] C.Kuo, C.Hanselman, Matlab tools for control system Analysis and design, New Jersey, Prentice Hall Internetional(UK), 1994, 223 p., ISBN 0-13-099946-6. [2] John Borrie, Modern control systems, New Jersey, Prentice Hall Internetional(UK), 1986, 317 p, ISBN 0-13-590290-8. [3] Katsuhiko Ogata, Modern control engineering, second edition, University of Minnesota, Prentice Hall Internetional,1990,951p., ISBN 0-13-598731-8 [4] Bruno Vanslambrouck, Trillingen en geluid, Hogeschool West-Vlaanderen, department PIH, 145 p. [5] Frederik D’hulster, Henk Lambrecht, Actieve trillingsonderdrukking bij hogesnelheidswerktuigmachines,Eindwerk, Faculteit Toegepaste Westenschappen, Departement Werktuigkunde, K.U.Leuven,1998, 97 p. [6] Mireille Deglorie, Identificatie van de mechanische aandrijving van een testbank, Eindwerk, departement PIH, Elektromechanica, 2003, 109 p. [7] Kurt Stockman, Analoge en digitale regeltechniek, Hogeschool West-Vlaanderen, department PIH. [8] System Identication Toolbox, Matlab User’s Guide, Matworks Inc., 1997. [9] Job Van Amerongen, Design in State Space,Lecture, Faculteit control engineering, universiteit Twente,2005.
81
Bijlage
1
Opmeten van de eigenfrequenties van de skilat.
bode 0
-10
-20
Amplitude
-30
-40
Reeks1
-50
-60
-70
-80 0
5
7,5
9
10
11
12
13
15
17,5
20
22,5
25
27,5
30
32,5
35
37,5
40
42,5
45
47,5
50
f(Hz)
2
Actuator V40
3
Versnellingsopnemer.
4
5
Commando’s matlab.
The Graphical User Interface Ident: Open the interface. Midprefs: Set directory where to store start-up information.
Simulation and Prediction Idinput: Generate input signals. Idsim: Simulate a general linear system. Pe: Compute prediction errors. poly2th: Create a model structure for input-output models defined as numerator and denominator polynomials. Predict: Compute predictions according to model.
Data Manipulation Dtrend: Remove trends from data. Idfilt: Filter data. Idresamp: Resample data.
Nonparametric Estimation Covf: Estimate covariance function. Cra: Estimate impulse response and covariance functions using correlation analysis. Etfe: Estimate spectra and transfer functions using direct Fourier techniques. Spa: Estimate spectra and transfer functions using spectral analysis.
Parameter Estimation Ar: Estimate AR model. Armax: Estimate ARMAX model. Arx: Estimate ARX model using least squares.
6
Bj: Estimate Box-Jenkins model. Canstart: Estimate multivariate models in canonical state-space form. Ivar: Estimate AR model using instrumental variable methods. Ivx: Estimate ARX model using general instruments. iv4: Estimate ARX model using four-stage instrumental variable method. Oe: Estimate Output-Error model. n4sid: Estimate state-space model using subspace method. Pem: Estimate general linear model.
Model Structure Creation arx2th: Define (multivariate) ARX structures. Canform: Generate canonical forms. mf2th: Create arbitrary linear model structure via an M-file that you write. Modstruc: Define state-space models with known and unknown parameters. ms2th: Create model structure for linear state-space models with known and unknown parameters. poly2th: Create a model structure for input-output models defined as numerator and denominator polynomials.
Manipulating Model Structures Fixpar: Fix parameters in structures to given values. Sett: Set the sampling interval. ss2th: Transform a state-space model to a parametrized canonical form. Thinit: Select or randomize initial parameter values. unfixpar Allow certain earlier fixed parameters be estimated.
7
Model Conversions Idmodred: Reduce a model to lower order. thc2thd :Transform from continuous to discrete time. thd2thc: Transform from discrete to continuous time. th2arx: Theta to ARX parameters. th2ff :Theta to frequency functions and spectra. th2par: Theta to estimated parameters and variances. th2poly: Theta to transfer function polynomials. th2ss :Theta to state-space matrices. th2tf :Theta to transfer functions. th2zp: Theta to zeros, poles, and static gains.
Model Presentation Bodeplot: Plot Bode diagrams. Ffplot: Plot frequency functions and spectra. Idplot: Display input-output data. Nyqplot: Plot Nyquist diagrams. Present: Display model on screen. Zpplot: Plot zeros and poles.
Information Extraction Getff: Extract the frequency functions from the freqfunc format. Get: Extract the sampling interval from the theta format. Getmfth: Extract the M-file name that defines the model structure. Getncap: Extract from the theta format the number of data upon which model is based. Getzp: Extract the zeros and poles from the zepo format. th2par: Extract estimated parameters and variances from the theta format.
8
Model Validation Compare: Compare model’s simulated or predicted output with actual output. Idsim: Simulate a model. Pe: Compute prediction errors. Predict: Predict future outputs. Resid: Compute and test model residuals.
Assessing Model Uncertainty Idsimsd: Simulate responses from several possible models. th2ff : Compute frequency function and its standard deviation.
th2zp: Compute zeros, poles, static gains, and their standard deviations.
Model Structure Selection Arxstruc: Compute loss functions for sets of ARX model structure. ivstruc: Compute loss functions for sets of output error model structures. Selstruc: Select structure. Struc: Generate sets of structures.
Recursive Parameter Estimation Rarmax: Estimate ARMAX or ARMA models recursively. rarx: Estimate ARX or AR models recursively. roe: Estimate Output-Error models (IIR-filters) recursively. Rpem: Estimate general input-output models using a recursive prediction error method. Segment: Segment data and estimate models for each segment.
9
1e orde bepaling ARX
u = meetwaarden.Y(1).Data' y= meetwaarden.Y(2).Data' een=iddata(y,u,0.005); een.InputName='insignaal'; een.OutputName='uitsignaal'; ze=een; plot(ze) ze=dtrend(ze); impulse(ze)
m2=arx(ze,[4 4 1]) compare(ze,m2) thred=idmodred(m2,1) ths=thd2thc(thred) A1 = ths.a B1 = ths.b C1 = ths.c D1 = ths.d [Ts,ns]=ss2tf(A1,B1,C1,D1) proces=tf(Ts,ns)
10
1e orde bepaling n4sid
u = meetwaarden.Y(1).Data' y= meetwaarden.Y(2).Data' een=iddata(y,u,0.005); een.InputName='insignaal'; een.OutputName='uitsignaal'; ze=een; plot(ze) ze=dtrend(ze); %impulse(ze)
m1=n4sid(ze,1,0.005) compare(ze,m1) thred=idmodred(m1) ths=thd2thc(thred) A1 = ths.a B1 = ths.b C1 = ths.c D1 = ths.d [Ts,ns]=ss2tf(A1,B1,C1,D1) proces=tf(Ts,ns)
11
2de orde ARX.
u = meetwaarden.Y(1).Data' y= meetwaarden.Y(2).Data' een=iddata(y,u,0.005); een.InputName='insignaal'; een.OutputName='uitsignaal'; ze=een; plot(ze) ze=dtrend(ze); impulse(ze)
m2=arx(ze,[12 12 2]) compare(ze,m2) thred=idmodred(m2,1) ths=thd2thc(thred) A1 = ths.a B1 = ths.b C1 = ths.c D1 = ths.d [Ts,ns]=ss2tf(A1,B1,C1,D1) proces=tf(Ts,ns)
12
2de orde n4sid.
u = meetwaarden.Y(1).Data' y= meetwaarden.Y(2).Data' een=iddata(y,u,0.005); een.InputName='insignaal'; een.OutputName='uitsignaal'; ze=een; plot(ze) ze=dtrend(ze); %impulse(ze)
m1=n4sid(ze,4,0.005) compare(ze,m1) thred=idmodred(m1) ths=thd2thc(thred) A1 = ths.a B1 = ths.b C1 = ths.c D1 = ths.d [Ts,ns]=ss2tf(A1,B1,C1,D1) proces=tf(Ts,ns)
13
Toestandsmodel.m: Hier worden de polen geplaatst op: -75.7±26i.
14
Generatebode.m: Ts=0.001 u1=u(500:9500); y1=y(500:9500); t1=t(500:9500); IN=frequentieanalyse(t1,u1,[0 1000 0 .5 -360 360],0); OUT=frequentieanalyse(t1,y1,[0 1000 0 .5 -360 360],0); N=length(IN)/2 fabs=(0:length(IN)-1)/(length(IN)*Ts); M_IN=abs(IN); P_IN=angle(IN); M_IN(1)=abs(IN(1)/2); M_OUT=abs(OUT) ; P_OUT=angle(OUT); M_OUT(1)=abs(OUT(1)/2); M=20*log10(M_OUT./M_IN); P=P_OUT-P_IN; w=fabs*2*pi; subplot(2,1,1),semilogx(w(2:2:length(w)),M(2:2:length(w))); title('BODE-DIAGRAM') xlabel('frequency(rad/s)') ylabel('magnitude') grid on subplot(2,1,2),semilogx(w(2:2:length(w)),rem(180/pi*unwrap( P(2:2:length(w))),360)); xlabel('frequency(rad/s)') ylabel('phase') grid on
15
Frequentieanalyse.m function X =frequentieanalyse (t,x,r,thresholdphase) fmin=r(1); fmax=r(2); amin=r(3); amax=r(4); pmin=r(5); pmax=r(6); n=length(t)-mod(length(t),2); X=fft(x,n); X=X*2/length(X); m=abs(X); p=angle(X); m(1)=abs(X(1)/2); N=length(X)/2; fabs=(0:length(X)-1):(2*N*0.01); subplot(3,1,1),var(fabs) AXIS([fmin fmax amin amax]) grid on xlabel('frequency [Hz]') ylabel('magnitude') title ('spectrumanalysis') for I=1:2*N, if abs(m(I))
16
17