Notat - Universitetet i Stavanger

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