Teoria sistemelor. Laborator 2: Liniarizarea ecuatiilor diferentiale.

Functii de transfer. Spatiul starilor. Raspunsul sistemelor.

Obiective: La sfarsitul acestui laborator studentii trebuie sa poata sa:

• Obtina o aproximare liniara a unei ecuatii diferentiale neliniare.• Simuleze un sistem descris prin ecuatii diferentiale.• Obtina functia de transfer pentru un sistem liniar.• Obtina un model ın spatiul starilor pentru un sistem liniar.• Calculeze si sa reprezinte grafic ın Matlab raspunsul unui sistem liniar la o intrare impuls sau treapta.

1 Liniarizarea ecuatiilor diferentiale

Exercitiul 1.1 Ecuatia Van der Pol

Se considera ecuatia Van der Pol neliniara:

d2x(t)

dt2+ (1− x(t)2)

dx(t)

dt+ x(t) = 0

saux(t) + (1− x(t)2)x(t) + x(t) = 0 (1)

a) Determinati o aproximare liniara a ecuatiei Van der Pol ın jurul punctului de echilibru x0 = 0, x0 =0, x0 = 0.b) Comparati solutia numerica a ecuatiei liniarizate cu cea a ecuatiei neliniare pentru un interval de timpt ∈ [0, 14]. Conditiile initiale sunt: x(0) = 0.001 si x(0) = 0.

Solutie:

a) Din ecuatia neliniara Van der Pol, pentru x0 = 0 si x0 = 0, se obtine x0 = 0.

Utilizati aproximarea prin serii Taylor pentru o functie de n variabile: y = g(x1, x2, ..., xn) ın jurul unuipunct x0 = (x10, x20, ..., xn0):

y = g(x10, x20, ..., xn0) +∂g

∂x1|x0

· (x1 − x10) + ...+∂g

∂xn|x0

· (xn − xn0)

In cazul unei ecuatii diferentiale, pentru seria Taylor considerati: x1 = x(t), x2 = x(t) and x3 = x(t) sifunctia neliniara g fiind membrul stang al ecuatiei neliniare Van der Pol equation. Se obtine:

g(x, x, x) = g(0, 0, 0) +∂g

∂x|(0,0,0)(x(t)− 0) +

∂g

∂x|(0,0,0)(x(t)− 0) +

∂g

∂x|(0,0,0)(x(t)− 0)

saug(x, x, x) = 0 + (−2xx+ 1)|(0,0,0)(x− 0) + (1− x2)|(0,0,0)(x− 0) + 1 · (x− 0)

g(x, x, x) = 0 + ∆x+∆x+∆x = 0

unde:∆x = x(t)− x0 = x(t)− 0

1

∆x = x(t)− x0 = x(t)− 0∆x = x(t)− x0 = x(t)− 0

Ecuatia:∆x+∆x+∆x = 0

saux(t) + x(t) + x(t) = 0 (2)

este liniara si omogena.

b) Vom determina acum solutia numerica a ecuatiei neliniare si a ecuatiei liniarizate, utilizand Matlabsi le vom compara. Considerati ecuatia neliniara Van der Pol (1). Vom transforma ecuatia diferentiala deordinul 2 ıntr-un sistem de doua ecuatii diferntiale de ordinul 1 (reprezentare ın spatiul starilor) si definimvariabilele de stare: x1(t) = x(t) and x2(t) = x(t) = x1(t). Ca urmare, x2(t) = x(t) si vom obtine modelulneliniar ın spatiul starilor:

x1(t) = x2(t)

x2(t) = −(1− x1(t)2)x2(t)− x1(t) (3)

Urmand aceeasi procedura, aproximarea liniara (2) va fi scrisa ca:

x1(t) = x2(t)

x2(t) = −x1(t)− x2(t) (4)

Conditiile initiale pentru ambele cazuri sunt: x1(0) = 0.001 si x2(0) = 0.

Matlab va calcula o solutie numerica a unui sistem de ecuatii diferentiale de ordinul 1 cu functiile ode23,ode45, ode15 sau altele (vedeti help pentru detalii). Pentru implementare ın Matlab este necesara o functiecare descrie membrul drept al ecuatiilor din sistem. De exemplu, sistemul (3) va fi descris ca o functieMatlab fnonlin.m care contine doar urmatoarele linii:

function f=fnonlin(t,x)

f=[x(2); -(1-x(1)^2)*x(2)-x(1)];

Sistemul (4) va fi scris ca o functie Matlab flin.m:

function f=flin(t,x)

f=[x(2); -x(1)-x(2)];

Programul principal care rezolva si compara solutia sistemului (3) si solutia sistemului (4) poate fi scrisca:

close all

clear all

x0=[0.001;0]; % initial conditions

[tn,xn]=ode23(@fnonlin,[0 14], x0);

[tl,xl]=ode23(@flin,[0 14], x0);

plot(tn, xn(:,1), ’b-’, tl, xl(:,1), ’r*’), grid on % solution x_1(t)

legend(’x_1 in VDP nonlinear’, ’x_1 in VDP linear’)

figure, plot(tn, xn(:,2), ’b-’, tl, xl(:,2), ’r*’), grid on % solution x_2(t)

legend(’x_2 in VDP nonlinear’, ’x_2 in VDP linear’)

Functia ode23 are cel putin trei argumente de intrare:

2

• @fnonlin: numele fuctiei care evalueaza membrul drept al ecuatiilor diferentiale

• [0 14]: un vector cu valorile timpului initial si final de integrare: [timp initial timp final]

• x0: un vector al conditiilor initiale

Functia va returna:

• tn: un vector coloana al valorilor de timp ıntre timpul initial si final

• xn: matricea solutiilor. Fiecare coloana din xn corespunde solutiei la timpul returnat in linia cores-punzatoare din tn.

Conditia initiala aleasa ın acest exemplu este foarte aproape de 0, astfel aproximarea ar trebui sa fiefoarte precisa.

Modificati conditiile initiale, pentru x1(0) = 0.5 sau x1(0) = 1 si comentati rezultatele.

Exercitiul 1.2 Pendulul simplu

Considerati pendulul mecanic din Figura 1.

Figura 1: Pendul

Pendulul consta dintr-un corp sferic mic si greu cu masa m suspendat de o tija rigida si foarte usoarade lungime l. Tija se poate roti ın jurul axei orizontale. Ca urmare, bila se roteste de-a lungul unui arcde cerc ın planul verical si pozitia sa este determinata de o singura coordonata, de exemplu, de deplasareaunghiulara notata cu x ın Figura 1. Miscarea bilei este guvernata de forta de greutate, mg, de forta defrecare Ff , si de momentul (cuplul) fortelor externe aplicate ın axa de rotatie , M(t).

Miscarea pendulului ın plan vertical este guvernata de ecuatia diferentiala:

ml2x(t) = M(t)−mgl sinx(t)− blx(t) (5)

unde:

• x(t) este unghiul pendulului fata de verticala• M(t) este momentul fortei ın punctul pivot• m este masa bilei, m = 0.5 kg• l este lungimea tijei, l = 1 m• g este acceleratia gravitationala, g = 9.8 m/s2

• b este coeficientul de frecare, b = 0.5

3

Considerati pendulul ca un sistem cu intrarea M(t) si iesirea x(t) si rezolvati urmatoarele probleme:

1. Determinati o aproximare liniara a ecuatiei (5) ın jurul x0 = 0, x0 = 0, x0 = 0 si pentru M0 = 0.

2. Utilizati aceeasi procedura ca ın aplicatia anterioara pentru a obtine si desena solutia numerica aecuatiei neliniare si a aproximarii ei liniare pentru M(t) = 0 si conditiile initiale:

(a) x(0) = π/4, x(0) = 0

(b) x(0) = π/2, x(0) = 0

Analizati si comentati rezultatele.

3. Determinati si reprezentati grafic solutia numerica a ecuatiilor neliniare si liniarizate, pentru conditiiinitiale nule si M(t) = 1.

4. Functia de transfer

(a) Din ecuatia diferentiala liniarizata, determinati functia de transfer de la intrarea M(t) la iesireax(t).

(b) Creati o functie de transfer ın Matlab utilizand functia tf si reprezentati grafic raspunsul siste-mului cu functia step.

(c) Utilizand functia Matlab impulse reprezentati grafic raspunsul sistemului la o intrare impuls.

5. Modelul ın spatiul starilor

(a) Din ecuatia liniarizata obtinuta la punctul 1 determinati un model ın spatiul starilor, cu variabilelede stare: x1(t) = x(t), x2(t) = x(t), intrarea u(t) = M(t) si iesirea y(t) = x1(t). Scrieti modelulın forma standard:

x(t) = Ax(t) +Bu(t)

y(t) = Cx(t) +Du(t)

(b) Creati modelul ın spatiul starilor ın Matlab utilizand functia ss si reprezentati grafic raspunsulsistemului la intrare impuls. Comparati rezultatul cu cel obtinut la punctul (4c).

Exercitiul 1.3 Trenuri Maglev

Trenurile Maglev functioneaza la viteza foarte mare si frecare foarte mica. Trenul leviteaza pe o distantade aproximativ 1 cm, dupa cum este aratat ın Figura 2. Trenurile construite pana acum utilizeaza douametode: prima utilizeaza forta de ”atractie” magnetica si a doua este suspensia electrodinamica utilizandbobine superconductoare pentru ”repulsie” magnetica.

Forta de levitatie FL este controlata de curentul i din bobinele de levitatie si poate fi aproximata prin:

FL = k ·i2(t)

z2(t)

unde z este distanta de levitatie. Aceasta forta este opusa fortei de greutate orientata ın jos F = mg.

Ecuatiile diferentiale care descriu miscarea pe verticala a trenurilor sunt scrise ca: suma fortelor = masa× acceleratia (vezi Figura 3)

a) mz(t) = mg − ki2(t)

z2(t)

b) mz(t) = ki2(t)

z2(t)−mg

4

Guidance magnets

Air gap z

Area of attraction

Air gap

z

Train car

Conducting plate

Stator coils

Rotor coil

Figura 2: Trenuri Maglev. a) Constructie cu atractie; b) Constructie cu repulsie

Figura 3: Fortele din sistem pentru cazul trenului: a) bazat pe atractie; b) bazat pe repulsie

1. Determinati relatia liniarizata ıntre pozitia verticala z si curentul i aproape de conditia de echilibru,pentru fiecare din cele doua cazuri. Pentru constante utilizati urmatoarele valori:

• Distanta de levitatie la echilibru z0 = 0.01 m• Masa trenului m = 10000 kg• Constanta fortei de levitatie k = 0.001 Nm2/A2

• Acceleratia gravitationala g = 10 m/s2

Idee: Trenul este la echilibru la distanta de levitatie nominala z0 = 1cm = 0.01 m. Astfel, pentruz0 = 0.01, acceleratia este zero: dz2/dt2 = z0 = 0, iar curentul de echilibru i0 poate fi determinat din:

mg − ki20z20

= 0.

Determinati o aproximare liniara pentru functiile:

f1(z, z, i) = mz(t)−mg + ki2(t)

z2(t)= 0

f2(z, z, i) = mz(t) +mg − ki2(t)

z2(t)= 0

ın jurul punctului (z0 = 0.01, z0 = 0, i0 = ...), ın termenii variatiilor: ∆z(t) = z(t) − z0, ∆i(t) =i(t)− i0, ∆z(t) = z(t)− z0.

2. Functia de transfer:

(a) Determinati functia de transfer de la curentul de intrare ∆i(t) la pozitia de iesire ∆z(t).Idee: Considerati sistemul cu intrarea ∆i(t) si iesirea ∆z(t) si utilizati ecuatia diferentiala lini-arizata pentru a calcula functia de transfer pentru fiecare din cele doua cazuri.

(b) Reprezentati grafic raspunsul la impuls al sistemelor ın bucla deschisa pentru un interval de timpde 0.1 secunde (cazul a)) si 1 secunda (cazul b)).Idee: Utilizati functiile Matlab impulse pentru a simula raspunsul la impuls al sistemului pentrufiecare caz ın parte. Analizati si explicati rezultatele.

3. Modelul ın spatiul starilor

5

(a) Din relatiile liniarizate obtinute la punctul 1, scrieti modelul ın spatiul starilor ın forma standard:

x(t) = Ax(t) +Bu(t)

y(t) = Cx(t) +Du(t)

Variabilele de stare sunt: x1(t) = ∆z(t), x2(t) = ∆z(t), intrarea ın sistem este curentul u(t) =∆i(t) si iesirea este pozitia verticala y(t) = ∆z(t).

(b) Creati un model ın spatiul starilor ın Matlab utilizand functia ss si reprezentati raspunsul laimpuls.

2 Raspunsul sistemelor si functii de transfer

Exercitiul 2.1

Un semnal de intrare r(t) = e−t este aplicat unui sistem cu functia de transfer G(s). Semnalul de iesirerezultat, pentru conditii initiale zero este:

c(t) = 2− 3e−t + 3e−2tcos2t, t ≥ 0

Determinati G(s) pentru acest sistem. (Utilizati Tabelul 1)

Exercitiul 2.2

Un sistem cu semnalul de intrare r(t) si de iesire y(t) este descris de ecuatia diferentiala liniara:

d2y(t)

dt2+ 2

dy(t)

dt+ 2y(t) = r(t)

1. Determinati functia de transfer, H(s).

2. Determinati polii si zerourile.

3. Determinati raspunsul sistemului la o intrare impuls, yi(t). (utilizati: y(t) = £−1{H(s)R(s)})

4. Determinati raspunsul sistemului la o intrare treapta unitara, ys(t).(Utilizati: y(t) = £−1{H(s)R(s)})

5. Utilizati functiile Matlab impulse si step pentru a desena raspunsul la impuls si treapta al sistemului.Comparati rezultatele cu yi(t) si ys(t), pentru un interval de timp t ∈ [0, 10].

f(t) L[f(t)] f(t) L[f(t)]

δ(t) 1 sin ata

s2 + a2

11

scos at

s

s2 + a2

t1

s2e−at sin bt

b

(s + a)2 + b2

e−at 1

s+ ae−at cos bt

s+ a

(s + a)2 + b2

Tabela 1: Tabel cu transformate Laplace

6