Reglarea temperaturii dintr-un cuptor cu rezistoare ventilate

Enunt: Se considera un sistem automat de reglare a temperaturii ca in figura de mai jos:

unde: reprezinta rezistenta termica care franeaza circulatia caldurii din interiorul cuptorului catre cavitatea de masurare;

reprezinta capacitatea calorica a cavitatii de masurare;

Figura 1

reprezinta temperaturile din interiorul si exteriorul cuptorului si din cavitatea de masurare;

reprezinta rezistenta termica care franeaza circulatia caldurii pornind de la surasa de caldura pana in interiorul cuptoului;

reprezinta capacitatea calorica a interiorului cuptorului; reprezinta rezistenta termica ce franeaza circulatia caldurii din interior spre exterior;

reprezinta capacitatea calorica exterioara presupunandu-se a fi infinita.

In acest caz temperatura este masurata prin intermediul unui termocuplu aflat intr-o teaca de masurare si a unui amplificator operational de instrumentatie care reda tensiunea . Caldura din cuptor vine de la rezistenta comandata in tensiune iar traductorul are caracteristica liniara.

Etape de proiectare:

I. Modelarea si modele

Se considera schema electrica a cuptorului cu rezistoare ventilate in figura 2:

In sistemul de mai sus , unde . In cazul schemei electrice caldura produsa devine:

Relatia de mai sus este notata cu (1). Expresia de legatura dintre temperatura

masurata si cea din interiorul cuptorului este: (2). Pe baza

celor doua relatii procesul este descris de urmatoarea schema functionala:

Q

aR

aC

mRm

mCeC

e

Figura 2

Modelul matematic va rezulta de forma:

Programul in Matlab pentru realizarea functiei de transfer si raspunsul sistemului la o intrare de tip treapta este urmatorul:

k1=100; k2=0.1;

ra=0.01; rm=3; rf=0.1;

ca=5000; cm=10;

tetam=rm*cm; tetafm=rf*cm; tetafa=rf*ca;

b0=k1*k2*rf/(rf*rm*cm*ca);

a0=1/(rf*rm*cm*ca);

a1=1/(rm*cm);

a=[0 1; -a0 -a1]; b=[0; b0]; c=[1 0]; d=0;

[num,den]=ss2tf(a,b,c,d,1);

sys=tf(num,den);

figure(1)

disp('Functia de transfer')

Figura 3

se

fR

1

1K+

+ sT1

sa sT2

sm2K

printsys(num,den);

step(sys)

grid;

num/den =

6.6667e-005

------------------------------

s^2 + 0.033333 s + 6.6667e-005

II. Proiectarea sistemului automat utilizand controlere PID

Pentru proiectarea regulatorului s-a realizat urmatoarea schema de reglare in bucla inchisa:

unde: - este functia de transfer pentru regulatorul clasic;

+

- sH R sH F

- este functia de transfer a partii fixate. Functia de transfer in bucla inchisa este:

(1)

In general functia de transfer in bucla inchisa de ordinul II este de forma:

(2)

Deoarece functia de transfer a partii fixate are forma (3)

atunci din relatiile (2) si (3) obtinem:

Regulatorul de tipP:

Se considera matricile rezultate din modelul matematic:

Se va utiliza o schema de reglare cu reactie dupa stare care are . In acest scop se verifica daca sistemul este complet controlabil,

adica se calculeaza rangul matricii m:

rangm=2 sistemul este complet

controlabil dupa stare si ecuatia caracteristica a sistemului este:

(ecuatia caracteristica reala)

Vom alege localizarile dorite ale polilor sistemului inchis. Dorim ca sistemul sa aiba un timp tranzitoriu ts sau tt rezonabil de mic si o amortizare rezonabila de la sistemele de ordinul 2 standard, alegand astfel polii doriti pentru sistemul inchis. La , unde

Observatie: Polii dominanti ai sistemului sunt si care au si pulsatia naturala . Ecuatia caracteristica dorita devine:

(ecuatia caracteristica dorita)

Pentru a determinamatricea k a amplificarilor cu reactie dupa stare se calculeaza cu relatia:

, unde matricea de transformare , unde m este matricea de controlabilitate, iar w se calculeaza cu ajutorul coeficientilor ecuatiei caracteristice reale astfel:

Matricea t devine:

iar

Matricea k a amplificarilor cu reactie dupa stare devine:

Semnalul de control u este determinat astfel:. Acest sistem este un sistem

regulator (de stare).

Schema in Simulink este urmatoarea:

%Proiectarea in Matlab:

a=[0 1; -a0 -a1]; b=[0; b0]; c=[1 0]; d=0;

a1=1/(rm*cm);

k1=100; k2=0.1;

ra=0.01; rm=3; rf=0.1;

ca=5000; cm=10;

tetam=rm*cm; tetafm=rf*cm; tetafa=rf*ca;

b0=k1*k2*rf/(rf*rm*cm*ca);

a0=1/(rf*rm*cm*ca);

[num,den]=ss2tf(a,b,c,d);

sys=tf(num,den);

m=[b a*b];

rang=rank(m);

a1=0.0333;

a2=0.0001;

alfa1=4;

alfa2=16;

w=[a1 1;1 0];

t=m*w; t1=inv(t);

k=[alfa2-a2 alfa1-a1]*t1;

In urma rularii s-a obtinut:

Proiectarea regulatorului de tip PI:

In acest caz se va folosi urmatoarea schema in Simulink:

Se considera matricile rezultate din modelul matematic:

Pentru regulatorul PI se vor calcula matricile:

si

semnalul de control: unde Se va determina amplificarilor cu reactie dupa stare prin utilizarea tehnicii plasarii polilor. Mai intai se examineaza rangul matricii P.

. Rangul acestei matrici este 3,

deci sistemul definit de ecuatia (*) este complet controlabil dupa stare si deci este posibila plasarea arbitrara a polilor. Ecuatia caracteristica pe sistemul dat de ecuatia (*) este:

(*)

Pentru a obtine performante rezonabile (viteza rezonabila si amortizare in raspunsul sistemului dorit) vom alege urmatoarii poli:

Ecuatia caracteristica corespunzatoare este:

In continuare se calculeaza matricea de transformare t data de relatia: unde

Inversa matricii t este:

Matricea va fi data de relatia:

%Proiectarea in Matlab

k1=100; k2=0.1;

ra=0.01; rm=3; rf=0.1;

ca=5000; cm=10;

tetam=rm*cm; tetafm=rf*cm; tetafa=rf*ca;

b0=k1*k2*rf/(rf*rm*cm*ca);

a0=1/(rf*rm*cm*ca);

a1=1/(rm*cm);

a=[0 1; -a0 -a1]; b=[0; b0]; c=[1 0]; d=0;

[num,den]=ss2tf(a,b,c,d);

sys=tf(num,den);

m=[b a*b a^2*b];

rang=rank(m);

a1=0.0333;

a2=0.0001;

a3=0;

alfa1=6;

alfa2=9;

alfa3=20;

w=[a2 a1 1;a1 1 0;1 0 0];

t=m*w; t1=inv(t);

k=[alfa3-a3 alfa2-a2 alfa1-a1]*t1;

III. Proiectarea sistemului automat cu ajutorul locului radacinilor

k1=100; k2=0.1;

ra=0.01; rm=3; rf=0.1;

ca=5000; cm=10;

tetam=rm*cm; tetafm=rf*cm; tetafa=rf*ca;

b0=k1*k2*rf/(rf*rm*cm*ca);

a0=1/(rf*rm*cm*ca);

a1=1/(rm*cm);

a=[0 1; -a0 -a1]; b=[0; b0]; c=[1 0]; d=0;

[num,den]=ss2tf(a,b,c,d);

sys=tf(num,den);

figure(1)

rlocus(sys);

IV. Determinarea caracteristicilor de frecventa Bode si a marginii de faza si a marginii de castig

k1=100; k2=0.1;

ra=0.01; rm=3; rf=0.1;

ca=5000; cm=10;

tetam=rm*cm; tetafm=rf*cm; tetafa=rf*ca;

b0=k1*k2*rf/(rf*rm*cm*ca);

a0=1/(rf*rm*cm*ca);

a1=1/(rm*cm);

a=[0 1; -a0 -a1]; b=[0; b0]; c=[1 0]; d=0;

%determinarea numitorului si numaratorului functiei de transfer

[num,den]=ss2tf(a,b,c,d);

%crearea vectorului de frecventa

w=logspace(-10,30,400);

%obtinerea amplitudinii si fazei pentru iesire

%trasarea caracteristicilor bode

%determinarea marginii de faza

figure(1)

%[mag,phase]=bode(num,den);

bode(num,den,w);

%margin(mag,phase,w);

grid;

%realizarea hodografului

figure(2)

%[re,im]=nyquist(num,den);

nyquist(num,den);

%plot(re,im);

title('Hodograful...');

Diagrama Bode si hodograful arata ca in figurile de mai jos:

V. Proiectarea in spatiul starilor a unui regulator LQ, LQR sau LQG

Termenul LQ este o prescurtare pentru regulatorul liniar optimal, cu indice de performanta patratic. Regulatorul LQ este cel mai important rezultat al controlului modern in abordare de stare.

La un regulator LQ performanta este estimata cu un indicator de performanta scalar de forma:

Unde R si Q sunt matrici simetrice pozitiv semidefinite sau pozitiv definite. Q>=0 si R>0.

In Matlab:

• K=lqr(A,B,Q,R) -rezolva problema regulatorului LQ liniar in timp continuu;

• Calculam matricea K a amplificarilor optimale de pe reactie de la stare astfel incat legea de control de la stare cu sa minimizeze J care este afectat ecuatiei cu constrangeri .

Trebuie sa se determine K=[k1 k2 k3 k4] astfel incat IP-ul sa fie minim. In acest caz q11 e suficient de mare comparativ cu q22, q33 si R alesi .

%Proiectarea in Matlab :

k1=100; k2=0.1;

ra=0.01; rm=3; rf=0.1;

ca=5000; cm=10;

tetam=rm*cm; tetafm=rf*cm; tetafa=rf*ca;

b0=k1*k2*rf/(rf*rm*cm*ca);

a0=1/(rf*rm*cm*ca);

a1=1/(rm*cm);

a=[0 1; -a0 -a1]; b=[0; b0]; c=[1 0]; d=0;

q=[100 0;0 1]; r=[0.01];

k=lqr(a,b,q,r);

k1=k(1); k2=k(2);

aa=a-b*k; bb=b*k1; cc=c; dd=d;

t=0:0.01:7;

%Obtinerea raspunsului sistemului automatizat LQ[y,x,t]=step(aa,bb,cc,dd,1,t);subplot(2,1,1)plot(t,y);grid;title('Raspunsul indicial al sist automatizat optimal LQ');xlabel('t sec');ylabel('iesirea y');

%Obtinerea curbelor x1, x2 functie de timpsubplot(2,1,2)plot(t,x);grid;title('Curbele x1, x2 functie de timp');xlabel('t sec');ylabel('x1, x2');

text(2.6,1.5,'x1');text(3,1,'x2');

S-a obtinut:

VI. Determinarea diagramei Nicols si a marginii de faza si a marginii de castig. Comparatie cu etapa IV

k1=100; k2=0.1;

ra=0.01; rm=3; rf=0.1;

ca=5000; cm=10;

tetam=rm*cm; tetafm=rf*cm; tetafa=rf*ca;

b0=k1*k2*rf/(rf*rm*cm*ca);

a0=1/(rf*rm*cm*ca);

a1=1/(rm*cm);

a=[0 1; -a0 -a1]; b=[0; b0]; c=[1 0]; d=0;

%determinarea numitorului si numaratorului functiei de transfer

[num,den]=ss2tf(a,b,c,d);

%crearea vectorului de frecventa

%w=logspace(-10,30,70);

%obtinerea amplitudinii si fazei pentru iesire

%trasarea caracteristicilor nichols

%determinarea marginii de faza

figure(1)

nichols(num,den);

margin(mag,phase,w);

Diagrama Nichols:

VII. Filtrul Kalman

k1=100; k2=0.1;

ra=0.01; rm=3; rf=0.1;

ca=5000; cm=10;

tetam=rm*cm; tetafm=rf*cm; tetafa=rf*ca;

b0=k1*k2*rf/(rf*rm*cm*ca);

a0=1/(rf*rm*cm*ca);

a1=1/(rm*cm);

a=[0 1; -a0 -a1]; b=[0; b0]; c=[1 0]; d=0;

te=1; gm=[1 0];

%discretizarea spatiului starilor

sys=ss(a,b,c,d);

sysd=c2d(sys,te, 'zoh');

[ad,bd,cd,dd,ts,td]=ssdata(sysd);

aa=[ad zeros(2,1);-cd*ad 1];

bb=[-bd;cd*bd];

cd=[1 0];

dd=0;

%estimator Kalman

q=0.1; r=0.1;

[kest,l,p]=kalman(sysd,q,r)

%comanda lqr

qq=eye(3);

rr=2000;

[k,s,l]=dlqr(aa,bb,qq,rr)

VIII. Proiectarea unui controler discret

Fie (1)

Ecuatia de comanda, care corespunde functiei de transfer in domeniul timp. In acest caz procedeul de discretizare utilizeaza metodele numerice uzuale pentru evaluarea componentelor integrale si derivative:

(2)

Daca utilizam metoda de extrapolare de ordin 0 avem:

(3)

Astfel, legea de comanda numerica pentru regulatorul PID este:

(4)

Este posibila utilizarea altor metode numerice pentru evaluarea componenetelor I si D, ca de exemplu:

(5)

Pentru metoda de extrapolare de ordin 1, si

(6)

In acest caz, legea de comanda numerica pentru regulatorul PID este:

(7)

Ecuatia (4) pentru momentul i-1 este:

(8)

Facand diferenta relatiilor (4) si (8) se obtine:

(9)

unde (10)Ecuatia (9) este forma incrementala a legii de comanda a regulatorului PID. Daca inlocuim si cu relatia (10) avem:

(11)

si (12)

unde (13)

Functia de transfer in cazul relatiei (12) este:

(14)

In cazul relatiei (7) functia de transfer este de forma:

(15)

In structura paralela a regulatorului, functia de transfer a componenatei I

numerice este: (16)

si daca utilizam expresia functiei de transfer a integratorului numeric, unde:

(17)

atunci folosim metoda Tustin de discretizare. Cu aceasta metoda, functia de transfer a componenetei derivative la limita

cauyala este: (18)

Unde (19) cu (20)

Schema regulatorului PID atunci cand comanda derivativa este data de functia de iesire este redata in figura 1, unde (21)

si (22)

Algoritmul de comanda al regulatorului PID este:

(23)

Se considera urmatoarea functie de transfer a regulatorului:

%Program in Matlab pentru discretizarea functiei de transfer a regulatorului

num =

4.2409 0.0317

zY c zE

zY

1

1

1

1

2

z

z

T

T

I

e

1

1

1

1

az

zb

+

+ -

zVI

zVD

RK zU

+-

>> den=[20.83 1]

den =

20.8300 1.0000

>> sc=tf(num,den) Transfer function:4.241 s + 0.03168----------------- 20.83 s + 1 >> sysd=c2d(sc,1,'zoh') Transfer function:0.2036 z - 0.2021----------------- z - 0.9531 Sampling time: 1

Functia de transfer discretizata a regulatorului este:

IX. Proiectarea unui controler robustX. Proiectarea unui controler fuzzy

lead=1;

lag=poly([-100 -400]);

sis=tf(lead,lag);

lsis=feedback(sys*sis,1);