Stavanger, 1. juli 2015 Det teknisknaturvitenskapelige fakultet ELE620 Systemidentifikasjon, 2015. Generell informasjon om faget er tilgjengelig fra fagets nettside, og for øvinger brukes It’s learning. Innhold 1 Tidsdiskrete signal 4 1.1 Frekvens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Nedfolding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Oppgave med tidsdiskret signal . . . . . . . . . . . . . . . . . . 7 2 Representasjon av system. 8 2.1 Differensligningen . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Differensialligningen . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 z-transferfunksjonen . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.1 Sammenheng med impulsresponsen . . . . . . . . . . . . 12 2.3.2 Stasjonær respons . . . . . . . . . . . . . . . . . . . . . . 12 2.3.3 Poler og nullpunkt . . . . . . . . . . . . . . . . . . . . . 13 2.4 s-transferfunksjonen . . . . . . . . . . . . . . . . . . . . . . . . 14 2.5 Diskret tilstandsrommodell . . . . . . . . . . . . . . . . . . . . . 15 2.6 Kontinuerlig tilstandsrommodell . . . . . . . . . . . . . . . . . . 16 2.6.1 Systemets egenverdier . . . . . . . . . . . . . . . . . . . 16 Karl Skretting, Institutt for data- og elektroteknikk (IDE), Universitetet i Stavanger (UiS), 4036 Stavanger. Sentralbord 51 83 10 00. Direkte 51 83 20 16. E-post: karl.skretting@uis.no. 3 Overganger 17 3.1 Linearisering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.2 Fra TRM til transferfunksjoner . . . . . . . . . . . . . . . . . . 18 3.3 Diskretisering av transferfunksjonen . . . . . . . . . . . . . . . . 19 3.4 3.3.1 Utledning av eksakt diskretisering . . . . . . . . . . . . . 20 3.3.2 Diskretisering med approksimasjoner . . . . . . . . . . . 21 3.3.3 Diskretisering med substitusjon . . . . . . . . . . . . . . 24 3.3.4 prewarp-metoden . . . . . . . . . . . . . . . . . . . . . . 26 3.3.5 Eksempel prewarp-metoden på sugefilter . . . . . . . . . 26 3.3.6 Poler ved diskretisering . . . . . . . . . . . . . . . . . . . 27 Diskretisering av TRM . . . . . . . . . . . . . . . . . . . . . . . 28 3.4.1 Eksakt diskretisering av lineær TRM . . . . . . . . . . . 28 3.4.2 Rekkeutvikling for diskretisering av lineær TRM . . . . . 30 3.4.3 Approksimasjon for diskretisering av TRM . . . . . . . . 31 3.4.4 Linearisering i kombinasjon med en diskretisering . . . . 32 4 Systemets egenskaper 4.1 4.2 Generell tidrespons . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.1.1 Delbrøkoppspalting for beregning av tidsrespons . . . . . 33 4.1.2 Matlab for beregning av tidsrespons . . . . . . . . . . . 34 4.1.3 Impuls- og steg(sprang)respons . . . . . . . . . . . . . . 34 Frekvensrespons . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.2.1 4.3 33 Utledningen av frekvensrespons . . . . . . . . . . . . . . 35 Stabilitet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 5 Noen viktige system 38 5.1 Forsterker . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5.2 Tidsforskjøvet signal . . . . . . . . . . . . . . . . . . . . . . . . 38 5.3 Integrator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2 5.3.1 Eksakt diskretisering av transferfunksjon . . . . . . . . . 40 5.3.2 Differensligning fra z-transferfunksjonen . . . . . . . . . 41 5.3.3 Differensligning ved differensialapproksimasjonen . . . . 41 5.3.4 z-transferfunksjonen fra differensligningen . . . . . . . . 41 5.4 Dobbelintegrator . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.5 Førsteordens system . . . . . . . . . . . . . . . . . . . . . . . . 43 5.6 5.5.1 Eksakt diskretisering . . . . . . . . . . . . . . . . . . . . 44 5.5.2 Eksempel med diskretisering i Matlab . . . . . . . . . . 45 5.5.3 Eksakt differensligning med tall . . . . . . . . . . . . . . 46 5.5.4 Diskretisering ved differensialapproksimasjon . . . . . . . 46 5.5.5 Eksempel med substitusjon . . . . . . . . . . . . . . . . 47 Andreordens system . . . . . . . . . . . . . . . . . . . . . . . . 48 5.6.1 Sprangrespons for andreordens system . . . . . . . . . . 50 6 Litt matematikk 53 6.1 Litt regning med komplekse tall . . . . . . . . . . . . . . . . . . 53 6.2 z-transformasjonen . . . . . . . . . . . . . . . . . . . . . . . . . 54 6.3 6.2.1 Eksempel med y(k) = ak . . . . . . . . . . . . . . . . . . 54 6.2.2 Eksempel med y(k) = kak . . . . . . . . . . . . . . . . . 55 6.2.3 Noen tranformasjonspar . . . . . . . . . . . . . . . . . . 55 6.2.4 Sluttverditeoremet . . . . . . . . . . . . . . . . . . . . . 56 Laplace-transformasjonen . . . . . . . . . . . . . . . . . . . . . 57 6.3.1 Eksempel . . . . . . . . . . . . . . . . . . . . . . . . . . 57 6.3.2 Noen tranformasjonspar . . . . . . . . . . . . . . . . . . 57 6.4 Om matriser . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6.5 Derivasjon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6.6 Om uttrykket eAt . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3 Notat om tidsdiskrete systemer. 1 Tidsdiskrete signal Her er litt repetisjon fra signalbehandlingen om frekvens, sampling og nedfolding. 1.1 Frekvens Frekvens (for et signal) er omdreininger, perioder eller gjentakelser per sekund. Enhet er Hertz, Hz. Frekvens kan også måles med vinkelhastighet, radianer per sekund, eller periodetid, antall sekund for en omdreining. Forholdet mellom disse enhetene er det viktig å ha klart for seg. Symbol Enhet Tekst Tp s f Hz ω rad/s Sammenheng periodetid T = 1 f = 2π ω frekvens f= 1 Tp = ω 2π vinkelhastighet ω = 2πf = 2π Tp Hvis det er klart ut fra sammenhengen at ω har enhet rad/s kan den ofte også kalles frekvens selv om ω egentlig er vinkelhastighet. Et signal med kun en frekvens, ω, er et sinussignal. For eksempel u(t) = sin(ωt + φ) = sin(2πf t + φ). (1.1) Her er φ en faseforskyvning og t er tid. 1.2 Sampling Sampling er måling av signalverdien ved regelmessige intervall, Ts (eller bare T ). Ts kalles tidssteg, sampleperiode eller sampletid. Som for signalet kan en også bruke samplefrekvens eller skjeldnere sample(vinkel)hastighet. Symbol Enhet Tekst Ts s fs Hz ωs rad/s Sammenheng tidssteg Ts = 1 fs = 2π ωs samplefrekvens fs = 1 Ts = ωs 2π samplehastighet ωs = 2πfs = 4 2π Ts En oppgir også i mange sammenhenger signalfrekvensen relativt til samplefrekvensen, for eksempel for frekvensrespons for digitale filter. Da bruker en ofte normalisert frekvens (=f /fs ) som er frekvens skalert slik at f = fs blir 1, eller som i Matlab (=2f /fs ) bruker en ofte at Nyquist-frekvensen (se del 1.3) fN blir 1. Et tredje alternativ er gjerne å skalere frekvensene (=2πf /fs ) slik at f = fs blir 2π og f = fN blir π. For eksempel med fs = 100 Hz, normalisert frekvens er da π/4 for f = 12.5 Hz. Når en bruker normalisert frekvens må en alltid oppgi hva fs (eller fN ) er normalisert til. 1.3 Nedfolding For å unngå nedfolding eller aliasing må signalet ikke inneholde frekvenskomponenter med frekvens over Nyquist-frekvensen, denne er halve samplefrekvensen: fN = fs /2 eller ωN = ωs /2. Nedfolding er at faktiske frekvenser høyere enn fN opptrer (speiles til) lavere frekvenser under samplingsprosessen. Ei Matlab-fil som illustrerer dette er fold.m. Det er viktig å vite hvordan nedfolding skjer. La f være signalets virkelige frekvens, signalet samples så med samplefrekvens fs . Da vil den tilsynelatende frekvensen etter sampling være f¯. Følgende oppskrift brukes for å finne f¯ gitt f og fs . 1. En finner først frekvensen f1 slik at f = f1 + mfs der m = 0, 1, 2, . . ., og 0 ≤ f1 < fs . Altså f1 = f − bf /fs cfs . 2. Hvis f1 ≤ fN = fs /2 så er f¯ = f1 , ellers, fN ≤ f1 ≤ fs , så er f¯ = fs − f1 . Her får en speiling. Et signal består ofte av flere frekvenskomponenter. Nedfolding vil gjelde alle frekvenskomponenter som har høyere frekvens enn Nyquist-frekvensen, fN = fs /2. Folding er illustrert i figur 1 og figur 2. Vi ønsker som regel ikke nedfolding (i forbindelse med regulering) og må derfor velge tilstrekkelig høy samplefrekvens, det vil si tilstrekkelig lite tidssteg. En løsning er å lavpassfiltrere signalene før sampling og på den måten forsikre oss om at signalene ikke inneholder frekvenskomponenter som kan foldes ned. Samplingsfrekvensen bør velges tilstrekkelig stor, en tommelfingerregel som er strengere enn Nyquist-kriteriet er ωs = 2π ≥ 10ω1 Ts der ω1 er høyeste frekvenskomponent i signalene, i.e. båndbredden. 5 (1.2) 6 1 4 ⇓ 1 2 3 4 5 4 1 3 2 7 4 6 -f /fs 2 normalisert frekvens, fs = 1 f¯ frekvens - fN fs Figur 1: Nedfolding, eksempel 1. Frekvenser mellom 0.6 og 0.9 med topp på 0.8 speiles ned til frekvenser mellom 0.1 og 0.4 med topp på 0.2. Frekvensen på 1.08 foldes ned til 0.08. 6 @ @ 1 4 ⇓ 1 2 3 4 5 4 1 6 fN 3 2 7 4 -f /fs 2 normalisert frekvens, fs = 1 f¯ frekvens - fs Figur 2: Nedfolding, eksempel 2. Frekvenser mellom 1.3 og 1.5 foldes til mellom 0.3 og 0.5. Frekvenser mellom 1.5 og 1.6 flyttes ned og speiles til å bli mellom 0.4 og 0.5. Frekvensen på 1.9 flyttes til 0.1. 6 1.4 Oppgave med tidsdiskret signal Gitt følgende signal y(t) = sin(4t + 1). a. Hva er frekvensen? b. Tidssteget, Ts , er 1 sekund. Hva er samplingsfrekvensen? c. Hvilken frekvens ses (etter sampling)? Svar: a. ω = 4 og f = ω 2π = 4 2π = 0.64 Hz. b. fs = 1/Ts = 1 Hz. c. fN = fs /2 = 0.5 Hz, og vi har her f > fN og får altså folding eller speiling. Obeservert frekvens blir f¯ = fs − f = 1 − 0.64 = 0.36 Hz. 7 u(k) - System h(z) y(k) - Figur 3: Enkel skisse av diskret LTI-system. Systemet er her gitt av ztransferfunksjonen h(z). Det er et inngangssignal (pådrag) u(k) og et utgangssignal y(k). 2 Representasjon av system. Et system er ofte en representasjon av en avgrenset fysisk prosess, eller en modell av en fysisk prosess. Generelt er system et noe videre begrep som ofte brukes om abstrakte, administrative eller logiske prosesser i videste forstand. Uansett blir systemet definert og avgrenset ut fra grensesnitt til omverda. Når systemet representerer en fysiske prosess er det signalene som danner grensesnittet i systemet (som tilsvarere krefter eller energioverføringer i fysisk prosess). Systemet er dermed relatert til signalene. Hvis signalene er tidsdiskrete så er også tilhørende system diskret, selv om det kan være en kontinuerlig fysisk prosess som ligger bak. Når (minst et av) signalene er kontinuerlige er også systemet kontinuerlig. Et system kan skisseres med en boks som i figur 3. Et blokkdiagram er en god måte å illustrere et system på. I et blokkdiagram brytes systemet opp i enklere elementer eller subsystem, disse tegnes ofte som blokker. Subsystemene forbindes med (interne) signaler, som tegnes som streker med retning (piler). Det aller enkleste eksempelet på et blokkdiagram er i figur 3, noen litt større system er i figur 10 side 44. Når systemet er bygd opp av enkle enheter, som en forstår godt, så kan en ofte bare ved å se på blokkdiagrammet forstå systemet og hva det gjør. Simulink, en tilleggspakke til Matlab, bruker blokkdiagram for å tegne systemet. Noen ulike måter å representere systemer på viser i figur 4. Representasjon av kontinuerlige system er til venstre og diskrete system til høyre. Systemet kan i tillegg representeres enten i tidsdomenet eller i transformert form, i henholdsvis s-domenet eller z-domenet. Når en modellerer prosessen ut fra de fysiske lover som gjelder så kommer en ofte fram til en eller flere differensialligninger, disse gir da en beskrivelse (modell) av systemet. Tilsvarende modellering av diskrete prosesser gir en beskrivelse (modell) med differensligninger. Differensialligningene kan ofte være ganske komplekse, og inkluderer både førstederiverte og høyere-ordens-deriverte og de er gjerne også ulineære. Ved å innføre nye interne signal, tilstander x, så kan alle differensialligninger uttrykkes med kun førstederiverte, og alle differensligninger med kun to etterfølgende 8 Kontinuerlig system - u(t) ẏ(t) = . . . Tidsdiskret system - - y(t) u(k) y(k + 1) = . . . 6 6 −1 L Z −1 Z L ? u(s) - - y(k) ? h(s) = y(s) u(s) y(s) - u(z) - h(z) = y(z) u(z) y(z) - Figur 4: Kontinuerlige og diskrete system. I de øverste boksene beskrives systemene med differensialligninger og differensligninger. I de nederste boksene beskrives systemene med transferfunksjoner. Laplace-transformasjonen og ztransformasjonen danner overgangen, per definisjon gjøres de på inngang- og utgangssignalene, men den enkleste formen fås oftest når transformasjonene gjøres direkte på ligningene og så ordnes. tidssteg, (k + 1) og (k). En tilstand er da en indre egenskap ved systemet, den kan være fysisk virkelighet slik som for eksempel en temperatur i systemet, men kan også være mer løsrevet fra en fysisk virkelighet. I første omgang er det gjerne enklest å tenke seg at x kan være nøyaktig det sammen som utgangen y, som når tilstand måles direkte. En beskrivelse med interne tilstander er en tilstandsrommodell. Fordelen med å innføre tilstander x er at en blir mer fleksibel i å definere systemet. Ikke minst fordi at ved å bruke hensiksmessige tilstandsvariabler så kan alle differensialligninger uttrykkes med kun førstederiverte. Det at en bruker flere tilstander kan også gi andre forenklinger i systembeskrivelsen. Ved å bruke tilstandsrommodell for differensialligninger (og differensligninger) så blir figur 4 slik som i figur 5. Her er tilstandsrommodellen i både generell og lineær form, og det er også med noe mer omkring systemene. Jeg synes denne figuren på en god måte viser de ulike måter å representere systemer på og den kan hjelpe å holde oversikten. Signalene, pådrag u, tilstand x og utgang eller måling y, kan gjerne være vektorer av flere signal. For eksempel for flere diskrete tilstander x1 (k) x2 (k) x(k) = x(k) = .. . . xn (k) Tilsvarende notasjon brukes også for u(k) og y(k), og for signalene i andre 9 domener, t, s og z. 2.1 Differensligningen Differensligningen viser hvordan utgangsverdi kan beregnes ut fra nåværende og tidligere innganger og tidligere utganger. Den generelle formen er y(k) = f (u(k), u(k − 1), . . . , y(k − 1), y(k − 2), . . .) (2.1) der f (·) er en funksjon (og har ingenting med frekvens å gjøre). Det er ofte mest hensiktsmessig å ha en lineær form på denne funksjonen. For et kausalt, lineært, tidsinvariant, diskret filter blir formen for f (·) en lineær differensligning. p r X X an y(k − n) (2.2) y(k) = bn u(k − n) − n=1 n=0 Her har en, slik det vanligvis gjøres, satt a0 = 1. 2.2 Differensialligningen Når en modellerer et system ut fra de fysiske lovene som gjelder så får en ofte en differensialligning, eller et sett med differensialligninger. Et enkelt eksempel er en integrator: ẏ(t) = Ki u(t). Konstanten Ki er en konstant. Det kan for eksempel være volum av væske i en tank, y(t) = V (t), som er integralet av strømmen inn, u(t) = qinn (t). Differensialligningen for et førsteordens system er 1 ẏ(t) = Ku(t) − y(t) T1 2.3 z-transferfunksjonen Boksen nede til høyre i figur 4 viser systemet representert med z-transferfunksjonen. Den er definert som forholdet mellom utgangssignalet og inngangssignalet sine z-transformasjoner h(z) = y(z) . u(z) (2.3) Hvis det er flere innganger og utganger så har systemet en z-transferfunksjon for hvert par av inngangs- og utgangssignal. h(z) vil da være en matrise av 10 Kontinuerlig Diskret Diskretisering −→ u(k) u(t)y(t) y(k) - h(z), x(k) h(s), x(t) u(s) y(s) u(z) y(z) Modell Fysisk system @ Matematisk modellering @ R @ D1 −→ ẋ = f (. . .) y = g(. . .) x(k + 1) = f (. . .) y(k) = g(. . .) Generell Tilstandsrommodell Lineær ↓ Linearisering ↓ D2 x(k + 1) = Φx(k) + Γu(k) ẋ = Ax + Bu −→ y = Dx + Eu ↓ L y(k) = Dx(k) + Eu(k) ↑ L−1 h(s) = y(s)/u(s) Z ↓ D3 −→ @ Observerte I @ @ Z −1 ↑ h(z) = y(z)/u(z) Transferfunksjon signal, u og y Eksperiment Figur 5: Oversikt over ulike systembeskrivelser og overganger. Det å forstå disse ulike måtene å beskrive systemer på, og ikke minst overganger mellom disse, er viktig i systemidentifikasjon. Her er signalet u pådraget, signalet y er (målt) utgang, systemets transferfunksjon er h, og systemets tilstand er x. En kan ofte forenkle med å måle tilstanden(e) direkte og da har en y = x. 11 transferfunksjoner, med dimensjon l × s der l er antall utganger og s er innganger. y1 (z) .. y1 (z) . · · · uy1s(z) u1 (z) (z) yl (z) y(z) y(z) .. .. . . = h(z) = = T = . . . u(z) u (z) u1 (z) · · · us (z) yl (z) yl (z) · · · us (z) u1 (z) antall z-transferfunksjonen representerer den direkte sammenheng mellom systemets inn- og utsignaler. z-transferfunksjonen kan finnes ut fra definisjonen over, eller ved å ta z-transformasjonen av differensligningen direkte. Dette siste gir ofte det enkleste uttrykket. En kan også finne z-transferfunksjonen ved diskretisering av s-transferfunksjonen. Merk: Når et kontinuerlig system diskretiseres blir signalene y(k) og u(k) avhengige av tidssteget T (eller Ts om en vil være helt presis i notasjonen). Dermed blir også y(z) og u(z) og z-transferfunksjonen h(z) avhengige av T . 2.3.1 Sammenheng med impulsresponsen Det er en tett sammenhengen mellom transferfunksjon og impulsrespons. Ligning 2.3 kan ordnes til y(z) = h(z)u(z). Hvis en lar pådraget være en enhetsimpuls, u(k) = δ(k) og u(z) = 1, gir det en alternativ definisjon: z-transferfunksjonen er z-transformasjonen av systemets impulsrespons, h(k) = yδ (k), h(z) = Z{h(k)} = yδ (z) 2.3.2 (2.4) Stasjonær respons Systemets stasjonære respons, også kalla stegrespons eller statisk respons ys , er den verdien utgangssignalet går mot når inngangssignalet er et sprang. Denne kan finnes ut fra sluttverditeoremet lim y(k) = ys = lim(z − 1)y(z) z→1 k→∞ (2.5) Sluttverditeoremet er enkelt å vise, se del 6.2.4. Den stasjonære responsen for et system h(z), det vil si et hvilket som helst system som er gitt av systemets transferfunksjon, er ut fra sluttverditeoremet ys = lim(z − 1)y(z) = lim(z − 1)h(z)u(z) z→1 z→1 Uz h(z) = lim U zh(z) = U hz (1). = lim(z − 1) z→1 z→1 z−1 12 (2.6) Jeg har her skrevet hz (1) for å få fram at det er h(z) som er innsatt 1, ikke impulsresponsen h(k). Når sprangets høyde ikke er definert kan en også sette U = 1. 2.3.3 Poler og nullpunkt Et (kausalt) lineært diskret system gitt med differensligningen p X an y(k − n) = n=0 r X bn u(k − n). n=0 z-transformene av signalene er y(z) = ∞ X y(k)z −k , u(z) = k=0 ∞ X u(k)z −k . k=0 En kan også ta z-transformasjon av koeffisientsekvensene a(z) = p X an z −n , b(z) = r X bn z −n . n=0 n=0 Transferfunksjonen for systemet er ∞ h(z) = y(z) b(z) X h(n)z −n . = = u(z) a(z) n=0 Transferfunksjonen for systemet kan alternativt skrives som y(z)a(z) = u(z)b(z), I begge tilfeller er a(z) og b(z) er polynom i z, eller z −1 om en vil. Vanligvis er gjerne polynomet i z −1 , det er da kausalt i forhold til differensligningen, da har en a(z) = a0 + a1 z −1 + a2 z −2 + · · · , og vanligvis har vi a0 = 1, og b(z) = b0 + b1 z −1 + b2 z −2 + · · · . z-transferfunksjonen for systemet er da h(z) = b(z)/a(z). Selv om polynomene oftes gis i z −1 , kan en også (ved å multiplisere både teller og nevner med en hensiktsmessig faktor z n ) uttrykke de som polynom i z. Dette kan være forvirrende, og det er viktig at en alltid har klart for seg hva polynomet er i, og rekkefølge for koeffisientene (hvilken koeffisient hører til hvilket ledd). Spesielt må en være oppmerksom på dette når en bruker Matlab, da dette ikke er gjort likt i alle funksjonene. For eksempel bruker filterfunksjonen A og B som polynomer i z −1 der første element er koeffisient for z 0 . Men dlsim-funksjonen1 har NUM og DEN som polynom i z med koeffisientene 1 dlsim funksjonen er i Control System Toolbox (CST) men er nå foreldet (obsolete), erstattet av lsim. 13 i avtagende rekkefølge, altså z 0 koeffeisienten for siste element. Kun når NUM og DEN har samme lengde virker denne funksjonen likt som filter. En representasjon av z-transferfunksjonen med poler og nullpunkt er ofte hensiktsmessig. Da faktoriseres hvert av polynomene a(z) og b(z) til faktorer der røttene inngår. Røttene i a-polynomet blir da polene og røttene i b-polynomet blir da nullpunktene for z-transferfunksjonen. Før faktorisering foretrekker jeg oftest å multiplisere med en hensiktsmessig faktor z n over og under brøkstreken, det vil si samme faktor for begge polynom, dermed er det polynom i z som faktoriseres. I Matlab kan en finne røttene for et polynom med roots-funksjonen. h(z) = (z − ζ1 )(z − ζ2 )(z − ζ3 ) · · · (z − ζr ) b(z) =K a(z) (z − z1 )(z − z2 )(z − z3 ) · · · (z − zn ) (2.7) • ζi er nullpunkt. • zi er poler. • K er en (forsteknings)faktor, K = b0 når a0 = 1. Generelt har vi at z-transferfunksjonen er utgang delt på inngang, (2.3). For en lineær differensligning kan z-transferfunksjonen ordnes slik at den skrives som en rasjonell funksjon, en brøk med to polynom, (2.7). y(z) og u(z) er ofte ikke polynom, men a(z) og b(z) er polynom. h(z) = b(z) y(z) = . u(z) a(z) (2.8) Systemets orden er den samme som grad for a(z) polynomet. 2.4 s-transferfunksjonen Boksen nede til venstre i figur 4 viser systemet representert med s-transferfunksjonen. Den er definert som forholdet mellom utgangssignalet og inngangssignalet sine Laplace-transformasjoner h(s) = y(s) . u(s) (2.9) En kan finne s-transferfunksjonen for et kontinuerlig system på tilsvarende måte som en finner z-transferfunksjonen for et tidsdiskret system, bare at en bruker Laplace-transformasjon i stedet for z-transformasjon. En kan altså ta Laplace-transformasjon av inngangs- og utgangssignaler, eller direkte av differensialligningene for systemet. 14 2.5 Diskret tilstandsrommodell Den diskrete tilstandsrommodellen er en generell ordnet form for differensligningene for systemet, tilsvarende er den kontinuerlige tilstandsrommodellen er en generell ordnet form for differensialligningene for systemet. Den viktigste utvidelsen er at en inkluderer systemets tilstander x. Dette kan ofte gi en mer fullstendig og korrekt systembeskrivelse. Med utgangspunkt i differensligninger får en da en diskret tilstandsrommodell, den er vist i boksen med x(k +1) i figur 5. Hver tilstand uttrykkes som en funksjon av tidligere tilstander og nåværende og tidligerer innganger. Videre kan hver utgang beskrives med en funksjon av nåværende og tidligere tilstander. En har ofte ikke med tidligere utganger eller nåværende og tidligere innganger ved beskrivelsen av utgangene, en kan ha det men en kan like gjerne bare ta med flere tilstander. En generell diskret tilstandsrommodell er på form x(k + 1) = f (x(k), u(k)) y(k) = g(x(k), u(k)) (2.10) Pådrag, måling, tilstand eller funksjonene kan gjerne være på vektorform. Med s pådrag, l målinger og n tilstander har en da u1 y1 x1 f1 g1 u2 y2 x2 f2 g2 u = .. , y = .. , x = .. , f = .. , g = .. . . . . . . us yl xn fn gl Den lineære tilstandsrommodellen kan skrives på noen ulike former, alt etter om en tar med alle ledd eller kun de viktigste. Her har vi et eksempel som inkluderer prosesstøy og målestøy x(k + 1) = Φx(k) + Γu(k) + Ωv(k), y(k) = Dx(k) + Eu(k) + w(k). (2.11) Signalene er som i den generelle diskrete tilstandsrommodellen. Matrisene er oppsummert i tabellen nedenfor Navn Symbol Dimensjon Transisjonsmatrise Φ n×n Pådragsmatrise Γ n×s Forstyrrelsematrise Ω n×n Målematrise D l×n Direktekoblingsmatrise E l×s 15 v(k) er stokastiske forstyrrelser på prosessen eller tilstandene, de går ofte gjennom matrisa Ω som ofte også er diagonal, men den trenger ikke være det. w(k) er stokastiske forstyrrelser direkte på hver enkelt måling, her har en ikke noen “blandingsmatrise” slik som Ω er for støyen v. 2.6 Kontinuerlig tilstandsrommodell For en kontinuerlig generell tilstandsrommodell, uten støyledd, så er den deriverte til hver tilstand en generell funksjon av alle tilstandene og alle pådragene. Tilsvarende er hver målingen en generell funksjon av alle tilstandene og alle pådragene. Dette kan skrives kompakt der alle symbol er vektorer ẋ = f (x, u) y = g(x, u) (2.12) Legg merke til at funksjonen f (·) her er forskjellig fra funksjonen f (·) i ligning 2.10. Her er f (·) en funksjon som gir den deriverte av tilstandene for en gitt tilstand og et gitt pådrag, mens i ligning 2.10 gir f (·) tilstanden i neste tidssteg. Det kommer som regel klart fram av sammenhengen hvilken funksjon en mener når en skriver f (·). En kontinuerlig lineær tilstandsrommodell, her uten støyledd, er ẋ = Ax + Bu y = Dx + Eu 2.6.1 (2.13) Systemets egenverdier Egenverdiene for det kontinuerlige systemet (2.13) er røttene i karakteristisk ligning det(sI − A) = 0. Røttene gir egenverdiene {si } = eig(A). Tilsvarende for diskret system (2.11) er egenverdiene røttene i karakteristisk ligning som her er det(zI − Φ) = 0. Røttene gir egenverdiene {zi } = eig(Φ). Ut fra sammenhengen Φ = eAT og egenverdidekomposisjon, A = XΛX −1 , så kan en vise at egenverdiene i det diskretiserte system henger sammen med egenverdiene i det tilsvarende kontinuerlige system {zi } = eig(Φ) = esi T (2.14) og omvendt {si } = eig(A) = 16 1 ln zi . T (2.15) 3 Overganger I systemidentifikasjon er det ofte viktig å kunne gå fra den en systemrepresentasjon til en annen. Det omhandler dette kapittelet. 3.1 Linearisering En (diskret) ikke-lineær prosess (system) kan ofte innenfor et begrenset operasjonsområde, i praksis omkring et arbeidspunkt, tilnærmes med god nøyaktighet med et lineært system. En har her gitt et system som er beskrevet ved en, eller et sett av flere, førsteordens differensligninger. Modellen er her som i (2.10) men vi innfører en forenklet notasjon for funksjonene x(k + 1) = f (x(k), u(k)) = f (k) y(k) = g(x(k), u(k)) = g(k) (3.1) Arbeidspunktet A er for tilstanden xA og pådraget uA . Utledningen her er med fast arbeidspunkt. Med varierende arbeidspunkt, A(k), blir utledningen ganske lik, den følger samme møster. Avvik fra arbeidspunktet skrives ∆x for tilstand og ∆u for pådrag. Tilstand og pådrag ved et tidspunkt er da x(k) = xA +∆x(k) og u(k) = uA + ∆u(k). Vi har x(k + 1) = xA + ∆x(k + 1) = f (k) = f [xA + ∆x(k), uA + ∆u(k)] Med Taylor rekkeutvikling av f (k) omkring arbeidspunktet får vi ∂f ∂f x(k + 1) = f (xA , uA ) + ∆x(k) + · · · + ∆u(k) + · · · ∂x A ∂u A Her har vi ikke skrevet ut leddene med høyere orden. Disse kan en anta blir små når en er nær arbeidspunktet, og dermed ser vi bort fra de. f (xA , uA ) blir jo tilstanden et steg etter arbeidspunktet, og en kan gjerne anta at en i løpet av et tidssteg fortsatt vil være ganske nær arbeidspunktet, altså f (xA , uA ) ≈ xA . Vi får da x(k + 1) = xA + ∆x(k + 1) ∂f ∂f xA + ∆x(k + 1) = xA + ∆x(k) + ∆u(k) ∂x A ∂u A ∂f ∂f ∆x(k + 1) = ∆x(k) + ∆u(k) ∂x A ∂u A ∆x(k + 1) = Φ∆x(k) + Γ∆u(k) ∂f ∂f Φ= Γ= ∂x A ∂u A 17 Merk at hvis arbeidspunkt er tidsvarierende så vil også matrisene Φ og Γ være tidsvarierende, altså Φ(k) og Γ(k). I ligningene over har en partiell derivert av en vektor med funksjoner med hensyn på en vektor med variabler. Dette trenger kanskje litt utdypning ∂f1 .. ∂f1 . ∂f1 · · · ∂x ∂x n 1 ∂fn ∂f ∂f . .. .. = Φ= = (3.2) = .. . . T ∂x ∂x ∂x1 · · · ∂xn ∂fn ∂fn · · · ∂xn ∂x1 Vektorene f og x er like lange, begge har lengde n systemet. Φ er da ei matrise med dimensjon n × n har en ∂f1 .. . ∂fn ∂f ∂f = Γ= = = ∂u ∂uT ∂u1 · · · ∂us som er antall tilstander i slik den skal være. For Γ ∂f1 ∂u1 .. . ∂fn ∂u1 ··· ... ∂f1 ∂us ··· ∂fn ∂us .. . (3.3) Γ blir da ei matrise med dimensjon n × s slik den skal være. Tilsvarende linearisering kan gjøres med måleligningen g(k) i (3.1) og en vil da få samme generelle uttrykk for D og E som for Φ og Γ bare en har g i stedet for f i (3.2) og (3.3). I praksis så er ofte D og E mye enklere enn Φ og Γ, ofte E = 0. 3.2 Fra TRM til transferfunksjoner Vi skal her gå fra diskret tilstandsrommodell, kapittel 2.5 (2.11), til transferfunksjon h(z). Tilstandsrommodellen er x(k + 1) = Φx(k) + Γu(k), y(k) = Dx(k) + Eu(k). Pådraget u, tilstanden x og målingen y er generelt vektorer. Direktekoblingsmatrisen E er ofte 0. Vi skal finne z-transferfunksjonen, generelt transfermatrisen h(z) med dimensjon l × s, fra u til y. Disse overgangene gjøres med transformasjoner. z-transformasjon av (2.11) gir zx(z) = Φx(z) + Γu(z) (zI − Φ)x(z) = Γu(z) x(z) = (zI − Φ)−1 Γu(z) 18 (3.4) og setter inn for x(z) i y(z) y(z) = Dx(z) + Eu(z) y(z) = D(zI − Φ)−1 Γu(z) + Eu(z) y(z) = D(zI − Φ)−1 Γ + E u(z) y(z) h(z) = = D(zI − Φ)−1 Γ + E u(z) Et eksempel der dette brukes er i del 5.4 side 42. 3.3 Diskretisering av transferfunksjonen Diskretisering vil si en overgang fra venstre del til høyre del i figurene 5 og 4. Vi starter med eksakt diskretisering, det gir et ikke helt enkelt uttrykk som kun er nyttig for enkle system. Eksakt diskretisering forutsetter et nullteordens sample- og holdeelement på inngangen. Approksimasjonene, i underkapitlene 2 og 3 med eksempel i 4, gir derimot noen enkle og nyttige substitusjonsregler. Til slutt her ser vi hva som skjer med polene ved diskretisering, de er viktige for systemets stabilitet. En kan også diskretisere tilstandsrommodellen, det skal vi se på i neste delkapittel. Et kontinerlig (fysisk) system kan diskretiseres med sampletida T , og da gi et tidsdiskret system, det vil si en diskret matematisk modell av systemet. Systemet er gjerne fysisk og konkret og har kontinuerlige egenskaper, for eksempel et lukka rom med en ovn der en måler temperaturen og pådraget (inngangen) er (elektrisk) effekt tilført ovnen. Når en diskretiserer systemet endrer en ikke på selve det fysiske systemet, men en endrer på signalene. For utgangssignalet så leses det nå kun av ved bestemte tidspunkt, nemlig ved tider t = kT . Men den fysiske egenskapen som måles, for eksempel temperatur, har selvsagt verdi også når den ikke måles. For inngangssignalet må en nå gjøre et valg for hvordan en vil at det skal virke. En må gjøre det diskrete signalet om til et kontinuerlig signal, dette gjøres med et holdeelement. Ofte kan en styre inngangssignalet, pådraget, og da er det rimelig å anta at en kun får endre verdier ved sampletidene, og en annen antakelse som også ofte gjøres er at den satte verdien da holdes konstant til neste gang verdien settes. Dette er et nullteordens sample- og holdeelement, det er gjerne en fysisk del av pådragsorganet. Selv når inngangssignalet er et kontinuerlig signal, som en ikke detaljert styrer, men kun måler ved sampletidene vil en ofte bruke et nullteordens sample- og holdeelement i modellen. En vil da bruke et tilnærmet signal inn til modellen, et trappesignal som er en tilnærmelse til det virkelige kontinuerlige signalet inn til systemet. Dette er ofte godt nok, og hvis det er for grov tilnærming kan en få det bedre med å ha kortere sampletid. Men det finnes alternativer 19 Kontinuerlig system h(s) - u(t) Tidsdiskret system =⇒ - h(z) - y(t) u(k) - y(k) ? Sa. h(z) u(k) ? δ(k) - Ho. - uh (t) h(s) y(t) - Sa. - y(k) Figur 6: Diskretisering av kontinuerlig system med sample- og holdeelement. som kan brukes hvis det ikke er råd å øke samplefrekvensen, det er Euler og Tustin (første og andre orden) og Runge-Kutta-metoden (høyere orden). I dette tilfellet er sample- og holdeelement en del av algoritmen i modellen. 3.3.1 Utledning av eksakt diskretisering Eksakt diskretisering er generelt vanskelig, men vi kan komme et stykke på vei, og det kan gjøres med (manuelle) beregninger i enkle tilfeller. Med velegnede verktøy, Matlab, kan eksakt diskretisering enkel gjøres for alle system. Vi har gitt systemets kontinuerlige transferfunksjon h(s), og skal finne h(z). Vi har jo at h(z) = y(z) når u(z) = 1, altså u(k) er en diskret enhetspuls. Starter med å finne kontinuerlig transferfunksjon for en diskret enhetspuls, 1 k=0 u(k) = δ(k) = (3.5) 0 ellers Når u(k) går gjennom et nullteordens sample- og holdeelement blir det kontinuerlige signalet, uh (t), som et rektangel. uh (t) = S(t) − S(t − T ). (3.6) der S(t) her er et enhetssprang. S(t) = 0 t<0 1 t≥0 Tar Laplace-transformasjonen av uh (t) (3.6) og får ved å bruke (6.30) uh (s) = 1 1 −sT 1 − e−sT − e = . s s s 20 (3.7) Merk at et hvilket som helst pådrag til systemet nå vil være et stykkevis konstant signal siden det er et nullteordens sample- og holdeelement på inngangen. Den kontinuerlige respons på en diskret enhetspuls for en vilkårlig system h(s) er 1 − e−sT y(s) = h(s)uh (s) = h(s) = (1 − e−sT )h(s)/s. (3.8) s Videre får vi tidsresponsen y(t) = L−1 {y(s)} = L−1 {(1 − e−sT )h(s)/s} (3.9) En må også ha verdi for yt (0) og vi lar denne være 0, yt (0) = 0. Sampling av y(t) gir y(k) og z-transformasjon av y(k) gir y(z) y(z) = Z{y(k)} = Z{yt (kT )} = ∞ X yt (kT )z −k k=0 = ∞ n X −1 −sT L {(1 − e )h(s)/s} |t=kT o z −k (3.10) k=0 Nå gjør vi et grep som kan virke litt merkelig, men som kan gjøres fordi alle operasjoner i formelen ovenfor er lineære. Faktoren e−sT i s-domenet er en forsinkelse med T i tidsdoment. Faktoren z −1 i z-domenet er en forsinkelse med ett tidssteg, som har lengde T . Faktoren (1 − e−sT ) kan da tas gjennom invers Laplace-transform til tidsdomenet, samples med tidssteg T til det diskrete domenet, og så tas z-transformasjon av til z-domenet der det nå blir til faktoren (1 − z −1 ). Merk også at h(z) = y(z) når u(z) = 1 som her. Altså får vi h(z) = y(z) = (1 − z −1 ) ∞ n o X L−1 {h(s)/s} |t=kT z −k , k=0 eller om en vil n o h(z) = y(z) = (1 − z −1 )Z L−1 {h(s)/s} |t=kT . (3.11) Lenger kommer vi ikke før vi har et uttrykk for h(s). Et eksempel der en bruker dette er for førsteordens system i kapittel 5.5.1. 3.3.2 Diskretisering med approksimasjoner Ved hjelp av differensial- og integralapproksimasjoner kan en forholdsvis enkelt diskretisere alle system gitt med (førsteordens) differensialligninger. En kan utlede enkle subtitusjonsregler som kan brukes på overgangen fra h(s) til h(z). En får ulike regler alt etter hvilken approksimasjon som brukes, men uansett er reglene enkle, de er i de etterfølgende delkapitlene, her tar vi kun selve approksimasjonene. 21 Figur 7: Differensialapproksimasjoner. Her illustreres at en kan tilnærme den deriverte for x(k) = xt (kT ) på flere måter. En kan gå forover ett tidssteg og bruke x(k) og x(k + 1) eller bakover og bruke x(k − 1) og x(k) eller en kan bruke sampleverdi både foran og bak. Stigningstallet for de rette linjer mellom samplingspunktene er tilnærminger til den deriverte. En fordel med metodene basert på approksimasjonene er at en gjerne kan ha et kontinuerlig inngangssignal u(t), ikke et stykkevis konstant signal slik som eksakt diskretisering krever. Vi bruker nå et tilstandssignal x(k) i stedet for utgangen y(k). Som vanlig bruker vi prikk-notasjon for den tidsderiverte ẋ = ẋ(t) = d d d x(t) og ẋ(k) = xt (kT ) = x(k) dt dt dt ẋ(t) kan gjerne være oppgitt som en funksjon av x(t) og u(t) og dermed blir ẋ(k) en funksjon av x(k) og u(k), denne funksjonen kan godt være ulineær også for lineære system. Men hvis en ikke har en eksplisitt formel for ẋ(k) så blir det vanskeligere siden ẋ(k) er en tidsderivert og kan ikke finnes ut fra kun de samplede punkta i x(k), en trenger x(t) for å finne eksakt verdi. Likevel kan ẋ(k) tilnærmes på flere måter ut fra nå-, forrige eller neste verdier. Dette er illustrert i figur 7. Tre nærliggende alternativer finnes, de gir Differensialapproksimasjonene som er 22 • Eulers forovermetode (EF) bruker x(k) og x(k + 1): ẋ(k) ≈ 1 x(k + 1) − x(k) . T (3.12) • Eulers bakovermetode (EB) bruker x(k − 1) og x(k): ẋ(k) ≈ 1 x(k) − x(k − 1) . T (3.13) • Tustins metode er gjennomsnitt av forover- og bakovermetoden. Tustins metode bruker x(k − 1) og x(k + 1). ẋ(k) ≈ 1 x(k + 1) − x(k − 1) . 2T (3.14) Integralapproksimasjoner brukes når en har ẋ(k), gjerne ut fra ẋ, gitt med et uttrykk eller en funksjon som en kan beregne. En kan da bruke dette til å finne eller approksimere x(k) når en har gitt x(k − 1). Utgangspunktet er ligningen Z kT x(k) = x(k − 1) + ẋ(t)dt (k−1)T Integralet her kan tilnærmes på tilsvarene måte som differensialene over. En har der satt ẋ(t) = f1 x(t), u(t) = f1 (t), en kunne like gjerne satt opp for eksempel Eulers forovermetode (EF) som x(k) ≈ x(k − 1) + T ẋ(k − 1). (3.15) Nå ser en tydelig at integralapproksimasjonen (3.15) egentlig bare er resultatet en får ved å snu litt på differensialapproksimasjonen (3.12). Tilsvarende kan en gjøre med EB (3.13) og får da x(k) ≈ x(k − 1) + T ẋ(k). (3.16) For Tustins metode kan en først legge merke til at en for differensialapproksimasjonen har at den er gjennomsnitt av Eulsers forover og Eulers bakover, altså (3.14) = 21 (3.12) + 12 (3.13). Vi vil beholde denne sammenhengen også for integralapproksimasjonen, altså (3.17) = 21 (3.15) + 12 (3.16), og dermed er x(k) ≈ x(k − 1) + T ẋ(k) + ẋ(k − 1) . 2 (3.17) Oftest brukes nok Eulers bakovermetode, fordelen er at den er enkel og at den kun bruker forrige og nåværende verdier for x(k) for å finne ẋ(k). 23 3.3.3 Diskretisering med substitusjon Ved hjelp av differensial- og integralapproksimasjoner kan en utlede enkle subtitusjonsregler som kan brukes på overgangen fra h(s) til h(z). En får ulike regler alt etter hvilken approksimasjon som brukes, men uansett er reglene gnaske enkle å bruke. La oss starte med Eulers forovermetode, der utgangspunktet er at den deriverte til et (kontinuerlig) signal ved et tidspunkt t = kT tilnærmes som i (3.12). Signalet kalles her x men det kunne like gjerne blitt kalla y. 1 x(k + 1) − x(k) . ẋ(t) = ẋ(kT ) ≈ T Substitusjonsreglene kan finnes vedå se på et eksempel med en enkel integrator. En integrator er gitt ved den kontinuerlige modellen (differensialligningen) ẋ = u Laplace-transformeres til sx(s) = u(s) og det gir transferfunksjonen h(s) = x(s)/u(s) = 1/s. Alternativ 1 Approksimasjon med Eulers forovermetode (3.15) gir x(k − 1) + T u(k − 1) T u(k − 1) T z −1 u(z) T z −1 u(z) T z −1 u(z) x(z) = 1 − z −1 x(z) T z −1 T h(z) = = = . −1 u(z) 1−z z−1 x(k) x(k) − x(k − 1) x(z) − z −1 x(z) x(z)(1 − z −1 ) = = = = (3.18) Alternativ 2 Approksimasjon med Eulers forovermetode (3.15) gir x(k + 1) x(k + 1) − x(k) zx(z) − x(z) x(z)(z − 1) = = = = x(k) + T u(k) T u(k) T u(z) T u(z) T x(z) = u(z) z−1 x(z) T h(z) = = . u(z) z−1 24 (3.19) For begge alternativer har en at h(s) = 1/s gir h(z) = T /(z − 1). Dette kan brukes til å gi følgende substitusjonsregel som kan brukes for en vilkårlig transferfunksjon h(s), men en må være klar over at dette er en tilnærmelse basert på Eulers forovermetode. h(z) = h(s)|s= z−1 (3.20) h(s) = h(z)|z=T s+1 (3.21) T Tilsvarende med Eulers bakovermetode (3.16) gir x(k) = x(k − 1) + T u(k) x(z) − z −1 x(z) = T u(z) T x(z) = u(z) 1 − z −1 T Tz x(z) = = . h(z) = −1 u(z) 1−z z−1 (3.22) Altså har en nå at h(s) = 1/s gir h(z) = T z/(z − 1) og får substitusjonsregler basert på Eulers bakovermetode h(z) = h(s)|s=(z−1)/T z h(s) = h(z)|z=1/(1−T s) (3.23) (3.24) Tustins metode (3.17) gir på samme måte T u(k) + u(k − 1) 2 T = u(z) + z −1 u(z) 2 T (1 + z −1 ) u(z) = 2(1 − z −1 ) T (1 + z −1 ) T (z + 1) = = . 2(1 − z −1 ) 2(z − 1) x(k) = x(k − 1) + x(z) − z −1 x(z) x(z) h(z) = x(z) u(z) (3.25) Dermed gir Tustins metode følgende substitusjonsregler h(z) = h(s)|s= 2(z−1) (3.26) h(s) = h(z)|z= 2+T s (3.27) T (z+1) 2−T s Substitusjonsreglene kan oppsummmeres i følgende tabell. EF h(s) → h(z) s= EB z−1 T s= h(z) → h(s) z = T s + 1 z = 25 z−1 Tz 1 1−T s Tustin s= z= 2(z−1) T (z+1) 2+T s 2−T s 3.3.4 prewarp-metoden Frekvenskorrigering er en modifisering av Tustins metode for diskretisering. Gitt et system ved s-transferfunksjonen h(s). Diskretiserer denne med Tustins approksimasjonsmetode med substitusjon, (3.26), og ser på frekvensresponsen etter diskretiseringen. Frekvensresponsen hz (ω) er z-transferfunksjonen innsatt z = ejωT . hz (ω) = hz (z)|z=ejωT = hs (s)|s= 2(z−1) |z=ejωT T (z+1) 2(ejωT − 1) = hs (s)|s= 2(ejωT −1) = hs T (ejωT + 1) T (ejωT +1) 2(ejωT − 1) · e−jωT /2 = hs T (ejωT + 1) · e−jωT /2 2 2j sin(ωT /2) 2j tan(ωT /2) = hs = hs T 2 cos(ωT /2) T 2 = hs (jv), der v = tan(ωT /2). T Vi husker at frekvensresponsen til et kontinuerlig system kan finnes ved å sette s = jv der v er frekvensen for innsignalet. Eksempelet i delkapittel 6.1 illustrerer nettopp uttrykket (z − 1)/(z + 1) i det komplekse plan. Det diskrete systemet hz (z) får altså samme frekvensrespons for frekvens ω som det kontinuerlige sysemet hs (s) får for frekvens v, der v = T2 tan(ωT /2). Merk at for små T så blir v ≈ ω. En kan dermed korrigere slik at det diskrete systemet, som jo er en tilnærming til det kontinuerlige systemet, gir eksakt samme respons som det kontinuerlige systemet ved en bestemt (kritisk) frekvens ω1 . Resultatet kan da oppsummeres med en (enkel) substitusjonsregel hz (z) = h(s)|s= 2ω1 (z−1) , T v1 (z+1) 3.3.5 v1 = 2 tan(ω1 T /2). T (3.28) Eksempel prewarp-metoden på sugefilter Her er et eksempel med et sugefilter, det er det samme som notch-filter. Vi tar her kun Matlab delen slik den blir med CST versjon 8.0. wf=1; T=1; zeta=sqrt(2)/2; teller = [1/(wf*wf), 0, 1]; nevner = [1/(wf*wf), 2*zeta/wf, 1]; sysc = tf(teller,nevner); % kontinuerlig system sysd = c2d(sysc,T,’tustin’); % diskret system 26 Figur 8: Bodediagram som viser frekvensresponsen for et sugefilter, det kontinuerlige tilfellet er med blå strek, diskretisering med Tustins metode er med grønn strek, og diskretisering med prewarp-metoden er med rød strek. sysdp = c2d(sysc,T,’prewarp’,wf); % diskret system med frekv.korr. bode(sysc,sysd,sysdp,linspace(0.1,10,500)); grid on; % Bode diagram for de tre systemene Resultatet viser i figur 8. 3.3.6 Poler ved diskretisering Poler for h(s) er s-verdier der nevneren er null, eller mer presist for en pol s1 har vi lim h(s) = ∞. s→s1 Diskretisering av h(s) gir forskjellige h(z) avhengig av metode, dermed får en også forskjellige poler avhengig av metode 27 Metode EF Substitusjon Pol z = 1 + Ts zi = 1 + T si EB z= 1 1−T s zi = 1 1−T si Tustin z= 2+T s 2−T s zi = 2+T si 2−T si zi = esi T Eksakt Substitusjonene i midtre kolonne er de en gjør i (3.20), (3.23) og (3.26). Eksakt diskretisering kan ikke gjøres med en substitusjon, men med nullte ordens holdeelement kom en fram til ligning 3.11 n o h(z) = (1 − z −1 )Z L−1 h(s)/s t=kT . og med regning som vi her ikke tar med, rekkeutvikling er gjerne enkleste måte, får en at en kontinuerlig pol si får tilsvarende diskrete pol zi . Sammenhengen mellom disse er gitt ved zi = esi T , si = 1 ln zi . T En merker seg at om det kontinuerlige systemet er ustabilt, poler i høyre halvplan, så blir også det (eksakte) diskretiserte systemet ustabilt med poler utenfor enhetssirkelen, altså Re(si ) > 0 ⇒ |zi | > 1. Substitusjonsreglene beholder ikke nødvendigvis stabiliteten! Et eksempel er med marginalt stabilt kontinerlig system som har pol på imaginær akse, for eksempel si = j. Eulers forovermetode gir da pol utenfor enhetssirkelen, ustabilt system. Eulers bakovermetode gir da pol innenfor enhetssirkelen, stabilt system. Mens Tustins metode gir her pol på enhetssirkelen slik som eksakt diskretisering også gir. 3.4 Diskretisering av TRM Som for diskretisering av transferfunksjonen kan også diskretisering av tilstandsrommodellen gjøres eksakt eller med en approksimasjon. Eksakt diskretisering gir også her et ikke helt enkelt uttrykk som kun er nyttig for enkle system. Diskretisering kan gjøres for differensialligningene, også når de er satt opp som tilstandsrommodell, enten den er generell eller lineær. 3.4.1 Eksakt diskretisering av lineær TRM Det vanligste er å diskretisere den lineære tilstandsrommodellen, det skal gjøres i dette kapittelet. Dette er en overgang fra systemet på form som i kapittel 2.6, (2.13). til form som i kapittel 2.5, (2.11). Lineær tilstandsrommodell, (2.13), 28 kan diskretiseres med nullteordens sample- og holdeelement, sampletid er T . Først skrives tilstandsligningen fra (2.13) som ẋ − Ax = Bu Med å multipliserer på begge sider med e−At får en e−At ẋ − e−At Ax ≡ d −At e x(t) = e−At Bu(t) dt Videre integrerer vi over et tidssteg fra t = t1 til t = t2 , (T = t2 − t1 ), noe som gir Z t2 −At t2 −At1 −At2 e−At Bu(t)dt x(t1 ) = x(t2 ) − e e x(t) t1 = e t1 Multipliserer så med e At2 og får A(t2 −t1 ) x(t2 ) − e t2 Z x(t1 ) = eA(t2 −t) Bu(t)dt t1 Her er integralet over u(t) i prinsippet varierende for t mellom t1 og t2 , med et nullteordens sample- og holdeelement så holdes u(t) = u(t1 ) i et tidssteg, altså helt fram til t2 . Dermed kan Bu(t1 ) trekkes utenfor integralet. Vi setter også x(t2 ) alene på venstre side og får Z t2 A(t2 −t1 ) x(t2 ) = e x(t1 ) + eA(t2 −t) dtBu(t1 ) t1 Med å sette t1 = kT , x(k) = xt (kT ) = x(t1 ) og t2 = (k + 1)T , x(k + 1) = x(t2 ) og skifte integrasjonsvariabel2 fra t til t2 − τ så blir dette (når en også bytter om integrasjonsgrensene) x(k + 1) = e AT Z T eAτ dτ Bu(k) x(k) + 0 Måleligningene blir uendret, det er ingen deriverte i den, og dermed har vi (2.11) x(k + 1) = Φx(k) + Γu(k), y(k) = Dx(k) + Eu(k). der Φ=e AT Z , og Γ = T eAτ dτ · B. (3.29) 0 2 Når en setter t → t2 − τ har en τ = t2 − t og dτ = −dt og integrasjonsgrensene blir t = t1 → τ = t2 − t = t2 − t1 = T og t = t2 → τ = t2 − t = 0 29 Alternativt kan Φ og Γ uttrykkes ved hjelp av invers Laplace-transform. Med å bruke ligning 6.42 får en n o Φ = eAT = L−1 L{eAt } t=T −1 −1 Φ = L (sI − A) (3.30) t=T For å finne tilsvarende R t uttrykk for 1Γ så bruker en en egenskap ved Laplacetransformasjonen, L{ 0 f (τ )dτ } = s L{f (t)}. Med å bruke dette så får en Z T nZ t o Aτ Γ = e dτ · B = eAτ dτ · B t=T 0 0 n nZ t o o eAτ dτ B = L−1 L t=T 0 n Z t o = L−1 L f (τ )dτ B med f (τ ) = eAτ t=T 0 o n n1 1 o = L−1 (sI − A)−1 B (3.31) Γ = L−1 L{eAτ } · B s s t=T t=T Matrisa B ble flyttet utenfor den innerste Laplace-transformasjonen. Det kan gjøres siden en antar at B er konstant, altså ikke en funksjon av tiden t, slik en også antar at A er konstant. 3.4.2 Rekkeutvikling for diskretisering av lineær TRM Formlene for Φ og Γ ovenfor er ikke helt enkle. Oftest, i praksis, tilnærmer en heller uttrykkene for Φ og Γ med rekkeutvikling, og en trenger sjelden mange ledd for å få ganske nøyaktig tilnærming. Har en mange ledd så kan en aksptere en stor T og likevel få god nøyaktighet. Men det kan likevel være nødvendig med liten T for mer kontroll med pådraget, eller bedre tilnærming til virkelig pådrag. 1. En regner først ut ei mellommatrise S = IT + 1 1 1 AT 2 + A2 T 3 + A3 T 4 + · · · 2! 3! 4! A er ei n × n matrise og I er ei n × n enhetsmatrise. T er her tidssteget, en skalar, den bør være tilstrekkelig liten, slik at de utelattet leddene ikke gir så stort bidrag. En tommelfingerregel er: T ≤ 2. Så regner en ut Φ med Φ = I + AS. 3. Og til slutt regner en ut Γ med Γ = SB. 30 0.1 |eig(A)|max (3.32) I Matlab, CST versjon 8.0, så er det c2d som brukes, denne funksjonen tar et LTI-system (ikke matriser for systemet) og metode som input. Et eksempel: A = [0,1; 0,0]; B = [0;1]; C = [1,0]; % C (og D) matrisene må også gis i ss Ts = 0.5; sysc = ss(A,B,C,0); % state space (tilstandsrommodell) sysd = c2d(sysc,Ts); % diskretiseres med tidssteg Ts 3.4.3 Approksimasjon for diskretisering av TRM Vi ser kun på tilfellet der tilstandsrommodellen approksimeres med Eulers forovermetode (EF). EF er enklest å utlede og ikke minst å bruke, sammenlignet med de vanligvis mer presise metodene Eulers bakovermetode og Tustins metode. Fra ligning 2.12 og ligning 2.13 har vi den kontinuerlige tilstandsrommodellen på henholdsvis generell og lineær form ẋ = f (x, u, . . .) ẋ = Ax + Bu + · · · (3.33) (3.34) Som vanlig brukes en notasjon der variabler uten indeks generelt er eller kan være vektorer eller matriser, og variabler eller funksjoner med indeks er skalarer. Dermed har en at (3.33) gjerne kan utvides til f1 (x(k), u(k), . . .) f1 (k) ẋ1 (k) .. .. ẋ(k) = ... = = . = f (k). . ẋn (k) fn (x(k), u(k), . . .) fn (k) Med Eulers forovermetode tilnærmes den deriverte som i (3.12) ẋ(k) ≈ 1 x(k + 1) − x(k) T Dette gir følgende diskretisering av modellen på lineær form (3.34), x(k + 1) = Φ x(k) + Γ u(k) x(k + 1) = x(k) + T ẋ(k) = x(k) + T Ax(k) + T Bu(k) = (I + AT )x(k) + BT u(k) (3.35) Merk at en i (3.35) bruker I, og ikke 1, siden x er vektor her. En ser også at dette tilsvarer rekkeutvikling når en kun har med til og med lineære ledd (S = IT , Φ = I + AS = I + AT og Γ = SB = BT ), det vil si utelater ledd med T 2 og høyere orden. 31 Modellen på generell form (3.33) kan gjerne diskretiseres direkte, x(k + 1) = x(k) + T ẋ(k) der en for ẋ(k) gjerne kan sette funksjonen f (x, u, . . .) fra ligning 3.33. 3.4.4 Linearisering i kombinasjon med en diskretisering En kan ta linearisering i kombinasjon med en diskretisering av TRM. En vanlig situasjon er at en ved å bruke matematisk modellering har fått systemet på form som en generell kontinuerlig tilstandsrommmodell, ligning 2.12. For å bruke Kalman-filter så ønsker en systemet på en form som en diskret lineære tilstandsrommmodell, ligning 2.11. Med Eulers forovermetode, ligning 3.12, får en at x(k + 1) = x(k) + T ẋ(k) ≡ f ∗ (k) Tidsteget er T . f ∗ (k) er her funksjonen som gir tilstand i neste tidssteg, som f (·) (uten *) i generell diskret modell, ligning 2.10 eller ligning 3.1. f ∗ (k) er brukt her for å skille tydelig fra funksjonen for ẋ i ligning 2.12 eller ligning 3.33. Vi kan nå finne matrisene Φ og Γ ut fra lineariseringen gitt i ligningene 3.2 og 3.3 med f ∗ (·) som funksjonen som gir tilstand i neste tidssteg. En får da, når en noen plasser ikke skriver (k), for lineariseringen etter Eulers forovermetode at ∂(x + T ẋ) ∂f ∗ Φ = Φ(k) = (3.36) = ∂x Ap ∂x Ap og ∂f ∗ ∂(x + T ẋ) Γ= (3.37) = ∂u Ap ∂u Ap Ap (merk ikke A som vi ofte ellers bruker for arbeidspunkt) er her arbedspunktet gitt av x og u. Ligningene over er tilsvarende som ligningene i kapittel 3.1. Merk at Eulers forovermetode kun brukes på de tidsderiverte og aldri på de partiellderiverte. 32 4 Systemets egenskaper 4.1 Generell tidrespons Hvis en har gitt systemet ved z-transferfunksjonen så kan en finne (simulere) tidsresponsen for et vilkårlig inngangssignal. Det er ikke nødvendigvis enkelt siden en først må finne z-transformasjonen av inngangssignalet og så den inverse z-transformasjonen av h(z)u(z). Dette kan en ofte løse med tabelloppslag for enkle funksjoner, med generelt er det vanskelig. 4.1.1 Delbrøkoppspalting for beregning av tidsrespons Delbrøkoppspalting må ofte gjøres for å finne tidsresponsen. Dette gjøres som regel i kombinasjon med tabelloppslag. Vi tar her et eksempel med delbrøksoppspalting der systemet er et førsteordens system h(z) = 0.3625z z − 0.8187 og inngangssignalet er et steg, med z-transform u(z) = 1 z−1 Utgangssignalet, i z-planet, er da y(z) = h(z)u(z) 0.3625z az bz = = + (z − 0.8187)(z − 1) z − 0.8187 z − 1 az(z − 1) + bz(z − 0.8187) az 2 − az + bz 2 − 0.8187bz = = (z − 0.8187)(z − 1) (z − 0.8187)(z − 1) 2 (a + b)z − (a + 0.8187b)z = (z − 0.8187)(z − 1) (4.1) som gir a + b = 0 og a + 0.8187b = −0.3625 med løsning a = −1.9994 og b = 1.9994, altså y(z) = z z 0.3625z = 1.9994( − ) (z − 0.8187)(z − 1) z − 1 z − 0.8187 Fra ligning 6.8 har vi at Z{ak } = z/(z − a), og dermed får vi her at y(k) = 1.9994(1k − 0.8187k ) ≈ 2(1 − 0.8187k ). 33 (4.2) 4.1.2 Matlab for beregning av tidsrespons Ofte er en enklere metode å finne ett og ett sample av utgangssignalet ut fra inngangssignalet, dette passer godt for datamaskiner. Dette kan gjøres med simulering i Matlab der slik simulering er gjort klar. En kan bruke funksjonene step og impulse og lsim når argumentet er et system. Mange klasser (objekttyper) kan brukes for å definere et system, i det vil si et lti-objekt3 . For eksempel: sys = tf(0.3625, [1, -0.8187], 0.2); step(sys); Tilsvarende som step har en impulse og den mer generelle lsim funksjonen. Merk at en også har definert funksjoner step og impulse for objekter av andre klasser (typer) enn LTI-objektene i CST, for eksempel iddata-klassen (og noen klasser til) i SIT. Generelt i Matlab må en vite hvordan argumentene til en funksjon brukes (tolkes) for å forstå resultatet riktig, det er spesielt viktig for disse funksjonene som kan ta polynom som argument. For eksempel når en bruker funksjonen sys = tf(teller,nevner,tidssteg) må en gi polynomkoeffisientene riktig, synkende grad i z. Det er alltid lurt å sjekke at det har blitt riktig, i Matlabkommandolinje skrives bare sys for å vise innhold i variabelen. 4.1.3 Impuls- og steg(sprang)respons Dette er systemets respons på noen enkle pådrag som ofte brukes for identifikasjon av systemer. Impulsresponsen er viktig på grunn av sammennhengen med systemets transferfunksjon. Stegresponsen er respons med pådrag som steg. For fysiske system er gjerne steg et relevant pådrag, mens impuls er mer problematisk. Statisk respons er det samme som stasjonær respons, stegresponsen etter lang tid. Eksempel med sprangrespons for andreordenssystem er i del 5.6.1. Hvis pådraget er statisk (konstant) så vil ofte også systemet få en statisk (konstant) tilstand. Med utgangspunkt i (2.11) uten støy v så kan den statiske respons finnes med å sette x(k + 1) = x(k) og u(k) = us , da må en ha xs = (I − Φ)−1 Γus . 4.2 (4.3) Frekvensrespons Utgangspunktet er et LTI system definert ved (impulsresponsen) h(z). Frekvensresponsen er det utgangssignalet en får når inngangssigalet er et sinussig3 Fra CST dokumentasjonen: An LTI object of the type TF, ZPK, SS, or FRD is created whenever you invoke the corresponding constructor function, tf, zpk, ss, or frd. 34 nal med en gitt frekvens (eller vinkelhastighet) ω. En kan vise at med inngangssignalet u(k) = sin(ωkT ) får en frekvensresponsen (utgangssignalet) y(z) = h(z)u(z). (4.4) Amplitude of fase(forskyvning) kan skrives som Amplitudefunksjonen er A(ω) = |hz (ejωT )| Fasefunksjonen er = ∠hz (ejωT ) Φ(ω) der • y(k) = A(ω) sin(ωkT + Φ(ω)) • |hz (ejωT )| er forsterkning eller dempning av amplituden, • ωkT er samme frekvens som inngangssignalet, og • ∠hz (ejωT ) gir faseforskyvning. 4.2.1 Utledningen av frekvensrespons Nå lar en ingangssignalet (pådraget) være et sinussignal med frekvens ω og amplitude U og fase φ = 0, sampletid er T . u(k) = U sin(ωkT ) = U Im(ejωkT ). (4.5) √ der j = −1. Vi skal nå finne utsignalet, det vil si frekvensresponsen. Det er litt regning og det krever at en er (godt) kjent med komplekse tall. y(k) = h(k) ∗ u(k) = ∞ X h(n)U sin[ω(k − n)T ] (konvolusjon) n=0 = U ∞ X jω(k−n)T h(n)Im[e ∞ X ] = U Im[ h(n)ejω(k−n)T ] n=0 n=0 ∞ X = U Im[ h(n)ejωkT e−jωnT ] y(k) = U Im[ n=0 ∞ X h(n)e−jωnT ejωkT ] (4.6) n=0 Merk at h(z) = P∞ n=0 h(n)z −n innsatt z = ejωT blir uttrykket i parentesen. hz (ejωT ) = ∞ X h(n)(ejωT )−n = n=0 ∞ X n=0 35 h(n)e−jωnT . Dermed får vi (4.6) kan skrives y(k) = U Im[hz (ejωT )ejωkT ] (4.7) Nå husker vi litt om multiplikasjon av komplekse tall. Sett at vi har to komplekse tall, z1 og z2 som multipliseres sammen z = z1 z2 . Vi har da at absoluttverdien av z blir |z| = |z1 ||z2 | og vinkelen eller fasen for z blir sum av vinklene for z1 og z2 , dette kan skrives ∠z = ∠z1 + ∠z2 . Skrives tallene med polare koordinater får en z1 = r1 ejθ1 , z2 = r2 ejθ1 , z = z1 z2 = (r1 r2 )ej(θ1 +θ2 ) = rejθ . (4.8) den imaginære delen av dette tallet, (z = z1 z2 ), er absoluttverdien multiplisert med sinus til vinkelen Im(z) = r sin θ = (r1 r2 ) sin(θ1 + θ2 ). (4.9) Det komplekse tallet i hakeparantesen i (4.7) er et produkt av to komplekse tall og resultatet ovenfor kan brukes. Husk at |ejωkT | = 1. Det gir jωT jωT y(k) = U hz (e ) · sin ωkT + ∠hz (e ) . (4.10) Merk at hz (ejωT ) er et komplekst tall uavhengig av k, gitt av transferfunksjonen, sampletida4 T , og frekvensen ω. Ut fra dette ser en at et sinussignal på inngangen gir et sinussignal på utgangen • |hz (ejωT )| er forsterkning eller dempning av amplituden • ωkT er samme frekvens som inngangssignalet • ∠hz (ejωT ) gir faseforskyvning Altså gir et lineært tidsinvariant (LTI) system en frekvensrespons som er en sinus med samme frekvens som inngangen, med med ei forsterkning eller dempning og med ei faseforskyvning som begge er gitt av transferfunksjonen h(z) innsatt z = ejωT , altså hz (ejωT ) der ω er inngangsfrekvensen. Amplitudefunksjonen er A(ω) = |hz (ejωT )| Fasefunksjonen er = ∠hz (ejωT ) Φ(ω) y(k) er rimeligvis avhengig av sampletida, og det er jo også u(k). Samplingsfrekvensen er fs = 1/T , en kan også velge å gi denne på samme måte som frekvensen for inngagssignalet, altså med en vinkelhastighet ω. Da har er ωs = 2π/T og tilhørende Nyquist-frekvens ωN = ωs /2 = π/T . Rimeligvis bør inngangssignalets frekvens være mindre enn Nyquist-frekvensen, altså ω < ωN = π/T eller T < π/ω. 4 Merk at hvis det diskrete systemet er funnet ved å diskretisere et kontinuerlig system så er også h(z) avhengig av sampletida T . 36 4.3 Stabilitet Stabilitetsegenskapen for et system kan defineres ut fra dets impulsrespons Asymptotisk stabilt system har at impulsresponsen dør ut, går mot null. For et diskret system har en limk→∞ yδ (k) = 0. En kan da vise at alle polene for h(z) er innenfor enhetssirkelen (ingen på). For en diskret tilstandsrommodell så må egenverdiene for matrisa Φ være innenfor enhetssirkelen. For et kontinuerlig system har en limy→∞ yδ (t) = 0. En kan da vise at alle polene for h(s) er i venstre halvplan. Marginalt stabilt system har at impulsresponsen fortsetter med samme størrelse, går mot en konstant (amplitude) ulik null men ikke uendelig. For et diskret system kan en da ha en eller flere (men ingen multiple) poler for h(z) på enhetssirkelen. For et kontinerlig system kan en da ha en eller flere (men ingen multiple) poler for h(s) på den imaginære aksen. Ustabilt system har at impulsresponsen går mot uendelig, i praksis krasjer mot ei fysisk grense. For et diskret system har en da minst en pol for h(z) utenfor enhetssirkelen eller multiple poler på enhetssirkelen. For et kontinuerlig system har en da minst en pol for h(s) i høyre halvplan eller multiple poler på den imaginære aksen. Matlab er velegnet for stabilitetsanalyse. 37 5 Noen viktige system Dette kapittelet analyserer noen enkle men viktige systemer. Vi kan her referere til figur 4 side 9 for å vise (de kontinuerlige og diskrete) signalene i forhold til (de kontinuerlige og diskrete) systemene. 5.1 Forsterker u(t) - K u(k) K I y(t)- - y(k) For en forsterker, blokkdiagram er skissert over, har en at den tegnes på samme måte i både diskret og kontinuerlig tilfelle. Forsterkning med en faktor K kan tegnes enten med en boks eller med en I. Transferfunksjonen er h(s) = K, h(z) = K (5.1) En helt enkel forstekning (eller dempning) av signalet. Ingen poler, ingen nullpunkt og ingen dynamikk. 5.2 Tidsforskjøvet signal u(t) - e−τ s y(t)- u(k) - z −n y(k)- Skissen ovenfor viser en forsinkelse med τ for det kontinuerlige systemet til venstre og en forsinkelse med n tidssteg for det diskrete systemet til høyre. Vi skal her kun se på det diskrete systemet. Merk at u(k) her er et vilkårlig inngangssignal med z-transformasjonen u(z). I det diskrete tidsdomenet gir et system med en tidsforsinkelse på n tidssteg y(k) = u(k − n). Negativ n her svarer til å forskyve signalet fram i tid. I praksis kan dette siste skje hvis tidsaksen erstattes med for eksempel en romlig akse eller hvis signalet er tatt opp på forhånd og avspilles ved bruk. 38 Vi skal nå finne z-transformasjonen av et tidsforskjøvet signal. y(z) = Z{y(k)} = Z{u(k − n)} ∞ X = u(k − n)z −k |m = k − n = k=0 ∞ X u(m)z −(m+n) m=−n = z −n ∞ X u(m)z −m (5.2) m=−n Merk at summen i (5.2) ikke er lik Z{u(k)} = u(z) siden summen ikke starter med leddet for m = 0. Hvis nå n ≥ 0, en forsinkelse, og vi antar (som vanlig) at u(m) = 0 for m < 0. Da får vi y(z) = Z{u(k − n)} = z −n ∞ X u(m)z −m = z −n u(z). (5.3) m=0 Hvis nå n < 0, altså en tidsforskyning framover, og vi har (som her) ensidig z-transformasjon må vi anta at de første (−n) sample av u(k) er null. Dette for å unngå at y(k), når samplene til u(k) skyves fram, ikke får verdier for negative indekser av k. Da får vi y(z) = Z{u(k − n)} = z −n ∞ X u(m)z −m = z −n u(z). (5.4) m=0 Altså, med de forutsetninger som er gjort, så får vi alltid (for alle n) Z{u(k − n)} = z −n Z{u(k)}. (5.5) z-transferfunksjonen for et system som forsinker inngangssiganel med n tidssteg, merk at hvis n < 0 så skyver systemet signalet framover, er da h(z) = y(z) = Z{u(k − n)}/Z{u(k)} = z −n . u(z) (5.6) z-transferfunksjonen for et diskret system med tidsforsinkelse på 1 eller n tidssteg er h(z) = z −1 , h(z) = z −n . (5.7) 39 Ki u(t) - Ki s y(t)- Q Q Q - - u(k) - T Ki z−1 T Ki + - z −1 I - y(k)- - 6 Figur 9: Blokkdiagram av en integrator med forsterkning Ki . Det kontinuerlige tilfellet øverst kan tegnes på to måter, symbolet brukt til høyre betyr det samme som boksen med 1/s inni. Forsterkning, med faktor Ki , tas ofte med sammen med integratoren. Det diskrete tilfellet nederst er en eksakt diskretisering, med nullteordens sample- og holdeelement, av det kontinuerlie systemet. 5.3 Integrator En integrator har den enkle differensialligningen ẏ = Ki u (5.8) Når en tar Laplace-transformasjonen av denne og ordner litt så får en stransferfunksjonen Ki y(s) = (5.9) h(s) = u(s) s 5.3.1 Eksakt diskretisering av transferfunksjon Eksakt diskretisering med nullteordens holdeelement gir med utgangspunkt i (3.11) n o h(s) h(z) = (1 − z −1 )Z L−1 { } |t=kT . s o n Ki h(z) = (1 − z −1 )Z L−1 { 2 } |t=kT . s ∞ X −1 −1 h(z) = Ki (1 − z )Z{kT } = T Ki (1 − z ) kz −k k=0 For å finne Z{k} brukes (6.9) med a = 1 som gir h(z) = T Ki (1 − z −1 ) z 1 = T Ki (1 − z −1 ) 2 (z − 1) (z − 1)(1 − z −1 ) 40 Dette gir T Ki T Ki z −1 = (5.10) z−1 1 − z −1 Merk at Eulers forovermetode her gir samme resultat mens EB og Tustin gir noen andre resultat. Felles for alle er likevel en pol i 1. h(z) = 5.3.2 Differensligning fra z-transferfunksjonen Ut fra z-transferfunksjonen h(z) = y(z) T Ki = z−1 u(z) får en (z − 1)y(z) = T Ki u(z) zy(z) = y(z) + T Ki u(z) ↓ Z −1 y(k + 1) = y(k) + T Ki u(k) (5.11) 5.3.3 Differensligning ved differensialapproksimasjonen Ovenfor har en brukt eksakt diskretisering og gått veien like til differensligningen. Ved å bruke differensialapproksimasjonen gitt av Eulers forovermetode (3.12) 1 ẏ(k) ≈ y(k + 1) − y(k) . T kan en gjøre dette mye enklere. En kan bruke dette direkte og får da ved å sette dette inn i differensialligningen (5.8) y(k + 1) = y(k) + T Ki u(k) (5.12) eller en kan bruke substitusjonsregel (3.20) i s-transferfunksjonen. Uansett, så er det basert på en approksimasjon, som for integratoren tilfeldigvis gir samme resultat som eksakt diskretisering. 5.3.4 z-transferfunksjonen fra differensligningen Har en gitt differensligningen (5.12) så kan en finne z-transferfunksjonen ved å ta z-transformasjonen av hvert enkelt ledd og så ordne resultatet på en 41 hensiktsmessig måte. y(k + 1) = y(k) + T Ki u(k) ↓ Z zy(z) = y(z) + T Ki u(z) (z − 1)y(z) = T Ki u(z) T Ki y(z) = u(z) = h(z)u(z) z−1 T Ki T Ki z −1 y(z) = = h(z) = . u(z) z−1 1 − z −1 5.4 (5.13) Dobbelintegrator En dobbelintegrator har den enkle differensialligningen ÿ = u (5.14) En kan sette denne opp som en lineær, kontinuerlig tilstandsrommodell ved å bruke to tilstander. En setter x1 = y og x2 = x˙1 = ẏ. Da har en x˙2 = ÿ = u. Tilstandsrommodellen er 0 ẋ1 0 1 ẋ = = x+ u ẋ2 0 0 1 y= 1 0 x Over er systemet på form som i kapittel 2.6, (2.13), der matrisene er 0 1 0 A= , B= , D= 1 0 , E = 0. 0 0 1 Vi skal nå diskretisere systemet slik som utledet i kapittel 3.4. Sampletida er T . Vi bruker formelen for den inverse av 2 × 2 matriser (6.31) der det trengst. Utgangspunktet for å finne Φ er (3.30) Φ = L−1 (sI − A)−1 t=T ( −1 ) s −1 = L−1 0 s t=T 1/s 1/s2 = L−1 0 1/s t=T 1 t 1 T = = (5.15) 0 1 t=T 0 1 42 Utgangspunktet for å finne Γ er (3.31) 1 o Γ = L (sI − A) B s t=T 1/s 1/s2 1 0 −1 = L 0 1/s s 1 t=T 2 3 1/s T /2 = = L−1 2 1/s T t=T −1 n −1 Dette gir altså den diskrete tilstandsrommodellen med matrisene 2 1 T T /2 Φ= , Γ= , D = [1, 0], E = 0. 0 1 T (5.16) (5.17) Transferfunksjonen h(z) kan nå finnes som utledet i kapittel 3.2. Vi finner da at z 0 1 T z − 1 −T zI − Φ = − = . 0 z 0 1 0 z−1 −1 1 z−1 T z − 1 −T −1 = (5.18) (zI − Φ) = 0 z−1 0 z−1 (z − 1)2 og h(z) = D(zI − Φ)−1 Γ + E 2 1 z−1 T T /2 +0 = [1, 0] 2 0 z−1 T (z − 1) 1 (z − 1)T 2 /2 + T 2 = [1, 0] (z − 1)T (z − 1)2 1 T 2 (z + 1) 2 2 2 = (zT /2 − T /2 + T ) = (z − 1)2 2(z − 1)2 T2 z −1 + z −2 = . 2 1 − 2z −1 + z −2 (5.19) Differensligningen er da lett å finne y(k) = 5.5 T 2 [u(k − 1) + u(k − 2)] + 2y(k − 1) − y(k − 2). 2 (5.20) Førsteordens system Differensialligningen for et førsteordens system skrives ofte som ẏ(t) = K 1 1 u(t) − y(t) = Ku(t) − y(t) . T1 T1 T1 43 (5.21) u u - - y- K T1 s+1 t(z) z−a y- u Q K +1/T1 ẏ Q Q + I I -6 u(n)- t(z) - + y(n + 1)- 6 z −1 y- y(n) - a J Figur 10: Blokkdiagram av et førsteordens system. Både blokkdiagram for kontinuerlig system (øverst) og blokkdiagram for diskret system (nederst) kan tegnes på flere måter. Når en tar Laplace-transformasjonen av dette og ordner litt får en s-transferfunksjonen K/T1 K = (5.22) h(s) = T1 s + 1 s + 1/T1 5.5.1 Eksakt diskretisering Vi har gitt transferfunksjonen for et førsteordens system h(s) = y(s) K K/T1 = = u(s) T1 s + 1 s + 1/T1 (5.23) Systemet skal diskretiseres. Det er nullteordens sample- og holdeelement på inngangen. Med å bruke (3.11) får en n h(z) = (1 − z −1 )Z L−1 o K/T1 . (s + 1/T1 )s t=kT (5.24) Vi setter c = 1/T1 og får videre, vi bruker (6.26), L−1 K/T1 c = KL−1 = K(1 − e−ct ) = K(1 − e−t/T1 ) (5.25) (s + 1/T1 )s (s + c)s 44 sampler med å sette inn t = kT n o o n K/T1 −t/T1 Z L−1 ) = Z K(1 − e t=kT (s + 1/T1 )s t=kT ∞ ∞ ∞ X X X = K (1 − e−kT /T1 )z −k = K z −k − e−kT /T1 z −k k=0 k=0 k=0 z z 1 1 = K − − = K 1 − z −1 1 − z −1 e−T /T1 z − 1 z − e−T /T1 Kz(1 − e−T /T1 ) z(z − e−T /T1 ) − z(z − 1) = . (5.26) = K (z − 1)(z − e−T /T1 ) (z − 1)(z − e−T /T1 ) Ut fra (5.24) og (5.26) får vi Kz(1 − e−T /T1 ) z − 1 Kz(1 − e−T /T1 ) = (z − 1)(z − e−T /T1 ) z (z − 1)(z − e−T /T1 ) 1 − e−T /T1 K(1 − e−T /T1 )z −1 = K (5.27) = z − e−T /T1 1 − e−T /T1 z −1 b bz −1 = . (5.28) h(z) = −1 1 − az z−a h(z) = (1 − z −1 ) Der en i (5.28) har satt b = K(1 − e−T /T1 ) = K(1 − a) og a = e−T /T1 . Eksakt diskretisering er nå gjort med å vise overgangen fra den kontinuerlige transferfunksjonen (5.22) til den diskrete transferfunksjonen (5.27). Orden til systemet er rimeligvis 1, det er det samme som antall poler (ulik 0) i systemet. For det diskrete systemet (5.28) er det en pol på den reelle aksen i 0 < a = e−T /T1 < 1 (siden T > 0 og T1 > 0), grensetilfellet med pol i 1 svarer til integrator, og grensetilfellet med pol i 0 svarer til en enkel forsinkelse. Telleren trenger ikke være en konstant for at systemet skal være et førsteordens system, en kan gjerne ha et tellerpolynom. I det diskrete tilfellet har en da generelt t(z) , 0 < a < 1. z−a Tellerpolynomet er ofte enkelt, og en skal være klar over at dette også i stor grad påvirker systemet. Tellerpolynomet svarer til et FIR-filter. Systemet kan implementeres i kaskade med utgangen av tellerpolynomet inn på et helt enkelt førsteordenssystem gjerne kun med en konstant over brøkstreken. h(z) = 5.5.2 Eksempel med diskretisering i Matlab Diskretisering med Matlab fungerer godt når en har tall for koeffisientene i teller og nevnerpolynomene for h(s). Vi har nå samme transferfunksjon og gitt verdiene K = 2, T1 = 3 [s] og T = 0.2 [s]. Det gir da a = e−T /T1 = e−1/15 ≈ 0.9355 og b = K(1 − e−T /T1 ) ≈ 0.1290. Dermed vil en eksakt diskretisering av h(s) = y(s) K 2 = = u(s) T1 s + 1 3s + 1 45 (5.29) med resultatet fra (5.28) være h(z) = b 0.1290 0.1290z −1 = = . z−a z − 0.9355 1 − 0.9355z −1 (5.30) I Matlab kan en finne dette med c2d-funksjonen i Control System Toolbox (CST). sysc = tf(2,[3,1]); sysd = c2d(sysc,0.2,’zoh’); og kan vise detaljene i objektet med get(sysd). Som alltid når en bruker polynomer i Matlab må en være helt klar over hvilke ledd de ulike koeffisienter hører til. 5.5.3 Eksakt differensligning med tall En kan finne den eksakte differensligningen ved invers z-transformasjon av z-transferfunksjonen (5.27). 0.1290z −1 y(z) = u(z) 1 − 0.9355z −1 y(z)(1 − 0.9355z −1 ) = u(z)0.1290z −1 ↓ Z −1 y(k) − 0.9355y(k − 1) = 0.1290u(k − 1) y(k) = 0.1290u(k − 1) + 0.9355y(k − 1) h(z) = 5.5.4 (5.31) (5.32) (5.33) (5.34) Diskretisering ved differensialapproksimasjon Ovenfor har en brukt eksakt diskretisering for å finne z-transferfunksjonen. Ved å bruke differensialapproksimasjonen gitt av Eulers forovermetode (3.12) kan en gjøre det samme. En kan starte med differensialligningen (5.21) direkte eller finne den med en invers Laplace-transform av s-transferfunksjonen (5.22) y(s) K = u(s) T1 s + 1 (T1 s + 1)y(s) = Ku(s) ↓ L−1 , T1 ẏ(t) + y(t) = Ku(t) h(s) = 46 (5.35) (5.36) yt (0) = 0 (5.37) Differensialligningen diskretiseres ved å sette t = kT , og å bruke differensialapproksimasjonen gitt av Eulers forovermetode (3.12) T1 y(k + 1) − y(k) + y(k) = Ku(k) T T1 T1 y(k + 1) + (1 − )y(k) = Ku(k) T T T T y(k + 1) + ( − 1)y(k) = Ku(k) T1 T1 T T y(k + 1) = Ku(k) + (1 − )y(k) T1 T1 (5.38) Dette gir ganske greitt differensligningen (5.38). Vi kan gjerne også ta z-transformasjon av dette og får T K T T T1 zy(z)+( −1)y(z) = Ku(z)y(z) = T1 T1 z − (1 − T ) T1 u(z) = 1− T Kz −1 T1 u(z) (1 − TT1 )z −1 og dermed z-transferfunksjonen h(z) = T Kz −1 y(z) T1 = . u(z) 1 − (1 − TT1 )z −1 (5.39) z-transferfunksjonen med EF (5.39) er nå ikke helt lik den eksakte z-transferfunksjonen (5.38). Bruker igjen verdiene som i Matlab eksempelet over, K = 2, T1 = 3 [s] og T = 0.2 [s]. Dermed blir differensligningen (5.39) y(k + 1) = 2 14 u(k) + y(k) = 0.1333u(k) + 0.9333y(k) 15 15 og z-transferfunksjonen (5.39) blir h(z) = 5.5.5 1 2 −1 z 15 − 14 z −1 15 = 0.1333z −1 . 1 − 0.9333z −1 Eksempel med substitusjon Vi tar et enkelt eksempel med substitusjon ut fra s-transferfunksjonen (5.22) og bruker Eulers forovermetode (3.20). Transferfunksjonen h(s) = y(s) K = u(s) T1 s + 1 47 (5.40) u - Kω02 2 0 +ω0 s2 +2ξω y- 2 Q u Kω0 ÿ - QQ + Q Q ẏ Q + I -6 + 6 y- 2ξω0 J ω02 J Figur 11: Blokkdiagram av et kontinuerlig andreordens system. skal nå diskretiseres med EF-substitusjon K z−1 T T1 s + 1 s= T T T K Kz −1 K T1 T1 = = T1 (z−1) z − 1 + TT1 1 − (1 − TT1 )z −1 +1 T h(z) = h(s)s= z−1 = = (5.41) Dette er jo ikke uventet det samme som (5.39) og når vi setter inn verdiene K = 2, T1 = 3 [s] og T = 0.2 [s] får vi selvsagt det samme resultat som i forrige del. 5.6 Andreordens system Andreordens system vil generelt gi svingninger, systemet har ofte en egenfrekvens som vil føre til forsterkning av et påtrykt sinussignal med samme frekvens. Det er derfor viktig med styring av inngangsignalet til et andreordenssystem. Tellerpolynomet vil selvsagt også være helt avgjørende for responsen, men en kan da dele systemet i to delsystemer i kaskade og dermed analysere de hver for seg. Sprangresponsen vil som regel gi innsvingninger mot en ny stasjonærverdi, eller hvis det er god demping, innsvingning uten oversving, indikasjonen på at det er andre og ikke førsteordenssystem er at andreordenssystem starter sprangresponsen mer forsiktig (flatere). s-transferfunksjonen skrives ofte på følgende form h(s) = Kω02 , s2 + 2ξω0 s + ω02 ω0 > 0, |ξ| < 1, (0 < ξ < 1). (5.42) ξ og ω0 og K er reelle tall. For at en skal få komplekskonjungerte poler i venstre halvplan (stabilt system) så må 0 < ξω0 og en velger da 0 < ξ og 0 < ω0 . Hvis |ξ| ≥ 1 så kan systemet skrives som en kaskade av to førsteordenssystem, og 48 s1 6 Im s J J ω J 0 J J J J α J J s2 s β φ, ξ = − cos φ Re - Figur 12: Skisse som viser sammenheng for ulike parametre som kan brukes for å beskrive et andreordens system. Parametrene inngår i ulike måte å angi de to komplekskonjungerte polene s1 og s2 på. analyseres enklest på den måten. Slik er det også hvis en i stedet for ω02 i nevneren skulle ha et negativt tall. En kan faktorisere nevnerpolynomet i transferfunksjonen h(s) = Kω02 p p (s + ω0 ξ − ω0 ξ 2 − 1)(s + ω0 ξ + ω0 ξ 2 − 1) (5.43) Kω02 (5.44) (s − s1 )(s − s2 ) der p s1 og s2 er ppolene i s-planet. Vi ser at disse er komplekskonjungerte. Siden ξ 2 − 1 = j 1 − ξ 2 kan polene nå skrives p p s1 = ω0 (−ξ + j 1 − ξ 2 ), s2 = ω0 (−ξ − j 1 − ξ 2 ). p med cos φ = −ξ, og da er sin φ = ± 1 − ξ 2 , kan vi skrive h(s) = s1 = ω0 ejφ , s2 = ω0 e−jφ . (5.45) Mer vanlig er det likevel å sette α = −ω0 ξ = ω0 cos φ, β = ω0 p 1 − ξ 2 = ω0 sin φ (5.46) Her ser en at α2 + β 2 = ω02 . Videre får en s1,2 = ω0 (cos φ ± j sin φ), 49 s1,2 = α ± jβ. (5.47) For å ha stabilt system med polene i venstre halvplan så må altså cos φ < 0 og dermed ξ > 0. En skisse som viser de to komplekskonjungerte polene i (venstre halvdel) av det komplekse s-planet er god for å vise sammenhengene mellom de ulike parametrene. ω0 er polenes avstand fra origo, α < 0 er reell del, β >p0 er imaginær del, og π/2 < φ < π er vinkel for pol, cos φ = −ξ og sin φ = 1 − ξ 2 . En fullstendig analyse av andreordenssystem tas ikke her. En kan bare vite at symbolene har følgende betydning (når polene er kompleks konjungerte par i venstre halvplan) • ω0 betegnes som udempet resonansfrekvens, det vil si den svingefrekvens en ville fått hvis dempingen α, eller ξ, var lik null. Merk at ω0 her ikke er det samme som knekkfrekvensen. • α er den absolutte dempningen. • ξ betegnes relativ dempingsfaktor. Den absolutte dempningen, α, angis da relativt til ω0 . ξ nær 0 er lite demping, sprangrespons vil da bli en svakt dempet sinuskurve, mens ξ nær 1 gir mye dempning (nesten som en førsteordens system). • φ er vinkel for pol og er mellom π/2 og π. Fra ligning 5.54 ser en at dette gir en faseforskyvning av sinusen i sprangresponsen. Med eksakt diskretisering og nullteordens holdeelement får vi i henhold til (3.3.6) de diskrete polene. z1,2 = es1,2 T = e(α±jβ)T = eαT e±jβT z1,2 = eαT cos(βT ) ± j sin(βT ) 5.6.1 (5.48) Sprangrespons for andreordens system Vi setter K = 1 i systemet i (5.44) og bruker α og β som i ligning 5.46, en har da ω02 = α2 + β 2 . Poler er som ligning 5.47. Sprang med størrelse 1 har Laplace-transformasjonen 1/s og vi får da at Laplace-transformasjonen for sprangresponsen blir y(s) = h(s)u(s) = α2 + β 2 h(s) = . s s (s − α − jβ) (s − α + jβ) (5.49) Her kan en ta delbrøkoppdeling og invers Laplace-transformasjon for å finne y(t). Det er noe arbeid, men ganske så rett fram, gjerne en god øvelse. Alternativt kan en se om en kan få hjelp av Symbolic Toolbox i Matlab, en må da gjerne prøve seg litt fram for å finne en lesbar og forståelig form. En skal 50 likevel være rimelig stø i manupulasjon av de matematiske funksjonene for å komme fram til et hensiktsmessig uttrykk. Jeg fant at følgende virket sånn noenlunde clear all j = sqrt(-1); syms s alpha beta s1 = alpha+j*beta; s2 = alpha-j*beta; ys = (alpha^2+beta^2)/(s*(s-s1)*(s-s2)); yt = ilaplace(ys); t1 = latex(ys) t2 = latex(simplify(yt)) Med R2014a fikk jeg da (med litt redigering) for t2: et (α−β j) (α2 + β 2 ) j et (α+β j) (α2 + β 2 ) j α2 + β 2 + − (α − β j) (α + β j) 2 β (α − β j) 2 β (α + β j) Med R2011a fikk jeg da (med litt redigering) for t2: − 1 − 2 β − j α et (α−β j) + j α et (α+β j) + β et (α−β j) + β et (α+β j) 2β og litt mer redigering gir 1 j α t (α+β j) 1− e − et (α−β j) − et (α−β j) + et (α+β j) 2β 2 og enda litt mer redigering gir jβt −jβt αt 1 jβt −jβt αt α 1 1 + e e −e − e e +e β 2j 2 = 1 + eαt α α sin(βt) − eαt cos(βt) = 1 − eαt cos(βt) − sin(βt) β β Resultatet for y(t) blir altså y(t) = 1 − eαt cos (βt) − α sin (βt) , β 0 ≤ t. (5.50) Første ledd er egentlig et enhetssprang, u(t), men her kan det settes til 1 siden u(t) = 1 for 0 ≤ t. Siden α < 0, har en at eαt går mot 0 når t vokser, og dermed er siste ledd en dempa sinus. Det at uttrykket inni parentesen er en (faseforskjøvet) sinus trenger kanskje litt forklaring. 51 En kan vise, kanskje enklest ved å uttrykke sinus og kosinusfunksjonene ved hjelp av eksponetialfunksjonen, at a sin x + b cos x = r sin(x + θ), r= √ a2 + b2 , tan θ = −b . a (5.51) Med a = −α/β og b = 1 og x = tβ blir (5.51) det samme som uttrykket i parentesen i (5.50). En får da r2 = 1 + r = p ξ2 1 ω02 ξ 2 = 1 + = 2 2 2 ω0 (1 − ξ ) 1−ξ 1 − ξ2 1 1 − ξ2 −b −1 β ω0 sin φ tan θ = = = = = tan φ a −α/β α ω0 cos φ θ = π−φ (5.52) (5.53) En må merke seg at når φ > π/2 (som her) så får en θ = π − φ og ikke bare θ = φ. Dermed y(t) = 1 − p eαt 1 − ξ2 sin βt + π − φ , 0 ≤ t. (5.54) Nå kan merke seg at yt (0) = 0 både for ligning 5.50 og ligning 5.54. En kan videre vise, kanskje enklest med utganspunkt i ligning 5.50 at y 0 (t) = ω02 αt e sin(βt), β 0 ≤ t. (5.55) Altså er y 0 (t) = 0 for t = kπ/β, k = 0, 1, 2, . . . og skifter også fortegn for disse t-verdiene. Men, spesielt hvis −α er (mye) større enn β så går eksponetialfaktoren så raskt mot null at den deriverte i praksis blir tilnærmet lik null før svingninger i y(t) kan observeres. Ut fra ligningene over kan en nå se hvordan sprangresponsen blir for ulike parametervalg. 52 Figur 13: Illustrasjon for regning med komplekse tall. 6 6.1 Litt matematikk Litt regning med komplekse tall Grunnleggende regning med komplekse tall forutsettes kjent. Det vil si representasjon med reell og imaginær del og med polar representasjon i det komplekse plan. Også de fire regnearter, polynomer og eksponesialfunksjonen med en del av egenskapene antas å være delvis kjent. Vi vil her bare forklare og illustrere noen egenskaper som brukes i dette dokumentet. Se på figur 13. Her her en markert et komplekst tall z på enhetssirkelen det vil si |z| = 1, z = ejωT = cos(ωT ) + j sin(ωT ). (6.1) Vi er nå interessert i absoluttverdien av z1 = z + 1. Vi legger merke til at de fire punkta: 0, z, z1 og 1 danner et parallellogram der alle sidene er like lange, |z1 | er lengden av den lengste (når ωT ≤ π/2) diagonalen i dette parallellogrammet, og at punktet der denne diagonalen skjærer enhetssirkelen dermed er z3 = ejωT /2 , merk at z3 ikke er midt på diagonalen. Multipliseres z1 med den konjungerte til z3 , z3∗ = e−jωT /2 , så får vi et reelt tall med samme absoluttverdi som z1 . |z + 1| = |z1 | = z1 z3∗ = (ejωT + 1)e−jωT /2 = ejωT /2 + e−jωT /2 = 2 cos(ωT /2). (6.2) 53 Ser vi på tallet z2 = z − 1 og multipliserer også dette med z3∗ = e−jωT /2 så får vi nå z2 z3∗ = (ejωT − 1)e−jωT /2 = ejωT /2 − e−jωT /2 = 2j sin(ωT /2). (6.3) Altså et rent imaginært tall. Dette viser at z2 og z3 står vinkelrett på hverandre, og at vinkel mellom z2 og imaginær akse er den samme som vinkel mellom z2 og reell akse, altså γ = ωT /2. Det komplekse tallet z4 er forholdet mellom z2 og z1 , vi har z4 = z2 z2 z3∗ 2j sin(ωT /2) = = = j tan(ωT /2) = j tan γ. ∗ z1 z1 z3 2 cos(ωT /2) (6.4) Altså med z = ejωT har en z−1 = j tan(ωT /2) = j tan γ. z+1 6.2 (6.5) z-transformasjonen Se gjerne Wikipedia sin artikkel, Z-transform, som også har en ganske så fullstendig tabell med transformasjonspar. For et diskret signal y(k) der k = 0, 1, 2, . . . kan en finne z-transformasjonen av signalet som ∞ X Z{y(k)} = y(z) = y(k)z −k . (6.6) k=0 Merk at vi har den ensidige z-transformasjonen her, det er oftest den en bruker i reguleringsteknikken da denne svarer til at en regner signalene lik null for negative indekser. 6.2.1 Eksempel med y(k) = ak Vi tar et nyttig eksempel her. Når a = 1 blir dette enhetssprang, u(k), og litt spesielt siden en har per definisjon at 00 = 1 så er dette impuls når a = 0. y(k) = ak , y(z) = Z{y(k)} = k = 0, 1, 2, . . . ∞ X k=0 ak z −k = (6.7) ∞ X a ( )k . z k=0 | az | Her må vi ha < 1 for å få konvergens, altsåP ROC for |z| > |a|. Da kan vi 1 k bruke formel for konvergent geometrisk rekke, ( ∞ |a| < 1). k=0 a = 1−a , y(z) = 1 1− a z = 1 z = . −1 1 − az z−a 54 (6.8) 6.2.2 Eksempel med y(k) = kak y(k) = kak , k = 0, 1, 2, . . . Fra forrige eksempel har vi Z{ak } = ∞ X ak z −k = k=0 z z−a Vi deriverer med hensyn på z, de to høgre uttrykkene må fortsatt være like hverandre ∞ X −a ak (−k)z −k−1 = (z − a)2 k=0 Ganger med −z på begge sider og får ∞ X k −k ka z k=0 az −1 az = = (z − a)2 (1 − az −1 )2 Vi ser at dette ut fra definisjonen av z-transformen gir følgende formel ∞ X Z kak = kak z −k = k=0 az −1 az = . (z − a)2 (1 − az −1 )2 (6.9) Her må vi ha ROC som i forrige eksempel siden det er det som er utgangspunktet og ingen nye restriksjoner er lagt til, altså ROC for |z| > |a|. 6.2.3 Noen tranformasjonspar Her er noen nyttige z-transformasjonspar, for en mer fullstendig tabell se Z-transform. En må være oppmerksom på konvergensområdet, Region Of Convergence eller ROC, når en finner z-transformasjonen. For Z{u(k)} bør en ha med at ROC : |z| > 1. Z{δ(k)} = 1, 1 z = Z{u(k)} = z−1 1 − z −1 δ(k) = 0 k 6= 0 1 k=0 der u(k) = 0 k<0 1 k≥0 (6.10) (6.11) Fra kapittel 5.2 fant vi z-transformasjon for en forsinkelse, h(z) = z −n . Bruker en der en enhetspuls som inngangssignal, altså δ(k) med Z{δ(k)} = 1, får en da at utgangssignalet blir δ(k − n) med z-transformen h(z)Z{δ(k)} Z{δ(k − n)} = z −n 55 (6.12) For ak og kak som ble utledet like foran her, så multipliserer vi med enhetsspranget for å presisere at sekvensene er 0 for k < 0, altså Z{ak u(k)} = 1 1− Z{kak u(k)} = 6.2.4 a z = 1 z . = −1 1 − az z−a az az −1 = . (z − a)2 (1 − az −1 )2 (6.13) (6.14) Sluttverditeoremet Gitt sekvensen x(k) der x(k) = xs , ∀ k ≥ M , der M < ∞ er et endelig heltall. Subscript s i xs står for stasjonær. Vi antar også at x(k) = 0, ∀ k < 0, slik vi gjør for ensidig z-transformasjon (kausale signal). z-transformasjonen av signalet er da ut fra definisjonen Z{x(k)} = x(z) = ∞ X x(k)z −k . (6.15) k=0 Merk at vi har ROC: |z| > 1 for kausal uendelig lang sekvens, her x(k). Vi definerer nå en ny sekvens: g(k) = x(k)−x(k −1), her får en da g(0) = x(0), g(1) = x(1) − x(0), g(2) = x(2) − x(1) og g(k) = 0, ∀ k > M . Siste element ulik null er da g(M ) = x(M ) − x(M − 1). g(k) har da endelig lengde og g(z) kan defineres for alle z, har ROC ∀z. En får da gz (1) = x(M ) = xs , altså stasjonærverdien for signalet x(k). Summen av de (k + 1) første leddene i sekvensen g(·) er da k X g(n) = x(0) + (x(1) − x(0)) + · · · + (x(k) − x(k − 1)) = x(k). (6.16) n=0 Og dermed får vi lim x(k) = lim k→∞ k→∞ k X n=0 g(n) = M X g(n) (6.17) n=0 siden g(k) = 0, ∀ k > M . For n = 0, 1, . . . , M er det ingen problemer med 1n = 11−n = 1 og heller ikke med limz→1 z n = limz→1 z 1−n = 1 og vi kan da sette, lim x(k) = k→∞ M X n=0 g(n) = M X M X 1−n g(n) · 1 = g(n) lim z n=0 z→1 n=0 (6.18) Bytter rekkefølge for grenseverdi og sum, og tar en faktor z utenfor summen og får M X lim x(k) = lim z g(n)z −n = lim z g(z) (6.19) k→∞ z→1 n=0 56 z→1 Siden g(k) = x(k) − x(k − 1) har vi z g(z) = z x(z) − z −1 x(z) = (z − 1) x(z) (6.20) merk at dette ikke gjelder for z = 1 men det gjelder for alle |z| > 1. Vi får da sluttverditeoremet med å kombinere de to siste ligningene lim x(k) = lim (z − 1) x(z). 6.3 (6.21) z→1 k→∞ Laplace-transformasjonen For mer enn definisjon og et eksempel av Laplace-transformasjonen se ei lærebok i matematikk eller Wikipedia sin artikkel Laplace transform. Der er det også mer fullstendige tabeller over transformasjonspar. Laplace-transformasjonen er definert som: Z ∞ L{f (t)} = F (s) = f (t)e−st dt (6.22) 0 6.3.1 Eksempel Vi tar nå et relevant eksempel på hvordan en kan finne Laplace-transformasjonen for en funksjon ut fra definisjonen Z ∞ −ct (1 − e−ct )e−st dt (6.23) L{1 − e } = 0 Z ∞ Z ∞ −st = 1e dt − e−ct e−st dt (6.24) 0 Z0 ∞ Z ∞ = e−st dt − e−(c+s)t dt 0 0 ↓ s > 0, s > −c h 1 i∞ h i∞ 1 −st −(c+s)t = − e − − e s c+s 0 0 1 c + s − s c 1 L{1 − e−ct } = − = = s c+s s(c + s) s(s + c) (6.25) (6.26) Ovenfor har vi egentlig funnet (utledet) to transformasjonspar, L{1} = 1s og 1 L{e−ct } = s+c . Laplace-transformasjonspar vil være oppgitt ved en eksamen, men de lineære egenskapene ved transformasjonen må ofte brukes hvis en skal transformere en bestemt funksjon. 6.3.2 Noen tranformasjonspar 1 L{1} = , s L{t} = 1 s2 og generelt L{tn−1 } = 57 (n − 1)! sn (6.27) L eat = 1 s−a L δ(t − a) = e−as e−as L u(t − a) = s 6.4 L teat = 1 (s − a)2 der δ(t) er en enhetsimpuls. 0 t<a der u(t − a) = 1 t≥a (6.28) (6.29) (6.30) Om matriser Ei 2 × 2 matrise og den inverse er 1 d −b a b −1 . A= , A = c d ad − bc −c a (6.31) Determinanten er: det A = ad − bc. Egenverdier for ei matrise er verdier λ slik at det(λI − A) = 0. Laplace-transformasjonen tas element for element i matriser. En har da for eksempel for ei 2 × 2 matrise med funksjoner, fi = fi (t) n f f o L{f } L{f } 1 2 1 2 L = (6.32) f3 f4 L{f3 } L{f4 } og for ei 2 × 2 matrise med funksjoner, hi = hi (s). n h h o L−1 {h } L−1 {h } 1 2 1 2 −1 L = h3 h4 L−1 {h3 } L−1 {h4 } (6.33) Tilsvarende, element for element regel, gjelder også for z-transformasjonen, for integraloperasjoner og differensialoperasjoner. 6.5 Derivasjon Partiell derivasjon av en vektor f med hensyn på en vektor av variabler x # " ∂f1 ∂f1 ∂f x1 f1 (·) ∂x2 1 = ∂x . (6.34) x= , f= gir ∂f2 ∂f2 x2 f2 (·) ∂x ∂x1 ∂x2 Trigonometriske funksjoner: d sin x = cos x dx d cos x = − sin x dx 58 (6.35) 6.6 Om uttrykket eAt Første gang en ser eksponentialfunksjonen med ei matrise som argument er gjerne litt rart. Egentlig er det ganske konsekvent, men det er gjerne nyttig å gå gjennom dette. Det bør være kjent at eksponentialfunksjonen ex og den komplekse varianten ez kan defineres med rekkeutvikling ∞ ex = 1 + x + X 1 1 2 1 3 x + x + ··· = xn . 2! 3! n! n=0 (6.36) En har x0 = 1 og følgende er definert for fakultet-funksjonen 0! = 1! = 1, n! = n(n − 1)! = 1 · 2 · 3 · · · (n − 1) · n (6.37) En kan definere eA , der A er ei kvadratisk matrise, tilsvarende som ex ∞ eA = I + A + X 1 1 2 1 3 A + A + ··· = An . 2! 3! n! n=0 (6.38) Her har en A0 = I der I er identitetsmatrisen med samme dimensjon som A. Ut fra ligning 6.38 ser en at eA også er ei matrise og den har sammen dimensjon som A. En vet jo at rekkeutviklingen av ex konvergerer for alle x, hva så med rekkeutviklingen av eA . Hvis A er ikke-defekt 5 , på engelsk “nondefective”, så har den en egenverdidekomposisjon A = XΛX −1 , det omvendte gjelder også slik at dette er en alternativ definisjon på ikke-defekt matrise. Λ er ei diagonal matrise med egenverdiene langs diagonalen, egenverdiene kan være komplekse også for reelle matriser. De fleste kvadratiske matriser har en egenverdidekomposisjon, selv om det også er uendelig mange matriser som ikke har det. Det er lett å vise at når A = XΛX −1 så konvergerer rekkeutviklingen A e eA ∞ ∞ ∞ X X 1 1 n X 1 −1 n A = (XΛX ) = XΛn X −1 (6.39) = n! n! n! n=0 n=0 n=0 n λ ∞ ∞ 1 X X 1 1 n −1 n −1 λ = X Λ X =X 2 X n! n! .. n=0 n=0 . e λ1 −1 eλ2 = XeΛ X −1 = X · (6.40) ·X ... Også når A er defekt så konvergerer rekkeutviklingen slik at eA er definert for alle kvadratiske matriser. 5 Ei matrise som har at alle egenverdier har geometrisk multiplisitet lik algebraisk multiplisitet er ikke-defekt, for mer fullstendig forklaring se ei lærebok i Lineær Algebra 59 Matlab-uttrykket e=exp(1); e^A beregner resultatet med formelen i ligning 6.40, den vil dermed ikke gi riktig resultat med defekte matriser. Matlab-funksjonen expm bruker en annen algoritme og finner eA også for defekte matriser, merk at Matlab-funksjonen exp med matrise-argument behandler hvert element i A for seg. Prøv for eksempel A=[2,1,0; 0,2,1; 0,0,2]; n=3; B=rand(n); e=exp(1); A1=e^A; A2=exp(A); A3=expm(A); B1=e^B; B2=exp(B); B3=expm(B); og se så på de resulterende matrisene. At er ei matrise der hvert element er en funksjon av t, med eksempel for 2 × 2 matrise, a11 a12 a11 t a12 t At = ·t= . (6.41) a21 a22 a21 t a22 t Laplace-transformasjonen for eAt er dermed kurant. Strengt tatt vises det her kun for ikke-defekte matriser, altså A = XΛX −1 . Vi bruker ligning 6.40, merk at da er At = (XΛX −1 )t = X(Λt)X −1 , og får L{eAt } = L{XeΛt X −1 } e λ1 t n eλ2 t = L X = X .. . −1 X L{eλ2 t } .. 1 s−λ1 . 1 s−λ2 60 o L{eλ1 t } = X .. −1 X . −1 X Videre kan vi nå ta den inverse av matrisene på begge sider og får 1 s−λ1 −1 1 −1 −1 At L{e } = X X s−λ2 .. . 1 −1 s−λ1 = X = X 1 s−λ2 ... X −1 s − λ1 s − λ2 .. . −1 X = X(sI − Λ)X −1 = XsIX −1 − XΛX −1 = sI − A Altså 1 (6.42) sI − A Det er samme form som det skalare uttrykket, som er når matrisa A har størrelse 1 × 1, 1 L{eat } = (s − a)−1 = (6.43) s−a L{eAt } = (sI − A)−1 = 61
© Copyright 2024