Politechnika Warszawska
Wydział Elektroniki i Technik Informacyjnych
Technika automatyzacji procesów
Projekt nr2
Wykonali:
inż. Adrian Tomasz Krężel
inż. Adrian Terwiński
inż. Piotr Szyszka
Prowadzący: dr inż. Sebastian Plamowski
Zasymulować działanie obiektu w Matlabie, a następnie:
Dobrać strukturę dwupętlowego układu regulacji z regulatorem PI lub PID bez odsprzęgania i z odsprzęganiem oraz nastawy tych regulatorów, dokonać analizy pracy zaproponowanego układu regulacji ,
Zaprojektować analityczny regulator predykcyjny, z uwzględnieniem ograniczeń przez rzutowanie, dokonać analizy pracy zaprojektowanego układu oraz porównać z dwupętlowym układem regulacji PI/PID,
Zaprojektować numeryczny regulator predykcyjny, z uwzględnieniem ograniczeń sterowania; Dokonać dokładnego porównania pracy układów regulacji z regulatorami predykcyjnymi numerycznym i analitycznym
W analizowanym układzie występują znaczne wzmocnienia obydwu wielkości sterujących na jedno wyjście regulowane.
Po analizie wartości wzmocnień zdecydowano się na regulowanie drugiej wartości regulowanej za pomocą pierwszego sterowania, oraz pierwszej zmiennej regulowanej za pomocą drugiego sterowania.
Po dołączeniu dwóch regulatorów do układu wyłączono człony inercyjne oraz całkujące obydwu regulatorów. Wzmocnienie członu proporcjonalnego jednego z regulatorów ustawiono na wartość 1, drugiego natomiast zwiększono aż do uzyskania niegasnącego przebiegu oscylacyjnego.
W ten sposób wyznaczono wzmocnienie krytyczne. Z wykresu odczytano okres oscylacji krytycznych i za pomocą współczynników Zieglera-Nicolsa wyznaczono K, Ti oraz Td regulatora PID.
Analogicznie wystrojono drugi regulator. Poniżej przedstawiono wykres niegasnących oscylacji drugiej wartości regulowanej oraz przebieg po wyznaczeniu wartości K, Ti oraz Td.
Po wyregulowaniu obiektu z dwoma regulatorami wykonano testy ze skokiem wartości zadanych w połowie symulacji. Początkowo wykonano skok wartości Ca z 260 [mol/m3] na 300 [mol/m3].
Następnie Ca ustawiono na wartości z punktu pracy i wykonano skok T z 393.9 [K] na 410 [K].
Konsekwencją dużych wzmocnień jednego sterowania na obydwie wartości regulowane jest duży peak, oraz oscylacje przy zmianie wartości zadanej którejkolwiek zmiennej.
W kolejnym kroku wyznaczono transmitancje odsprzęgajace wpływ działania pierwszego sterowania na drugie wejście oraz drugiego sterowania na pierwsze.
Stosując odsprzęganie otrzymano następujące przebiegi wartości regulowanych:
Na otrzymanym obiekcie wykonano analogiczne testy co na obiekcie bez odsprzęgania. Skok wartości Ca:
Skok wartości zadanej T:
Po zastosowaniu odsprzęgania uzyskano mniejsze przeregulowania obiektu. Przebiegi wartości regulowanych są też mniej oscylacyjne, bardziej stabilne.
Schemat badanego układu składającego się z obiektu oraz regulatora MPCS
Obiekt regulacji utworzono na podstawie zlinearyzowanych macierzy stanu A, B , C modelu, który obliczono w pierwszej części projektu.
$$A = \begin{bmatrix}
- 7,5410 & - 91,3059 \\
0,8503 & 5,5280 \\
\end{bmatrix}$$
$$B = \begin{bmatrix}
1 & 0 \\
0 & - 6,0620 \\
\end{bmatrix}$$
$$C = \begin{bmatrix}
1 & 0 \\
0 & 1 \\
\end{bmatrix}$$
Prawo regulacji regulatora MPCS wyznaczono korzystając ze wzoru:
x(k) – stan modelu z aktualnej chwili czasowej
u(k-1) – sterowanie z poprzedniej chwili czasowej
Yzad(k) – wartość zadana wyjść regulowanych
Macierz ${\overset{\overline{}}{K}}_{1}$ stosowana we wzorze na prawo regulacji regulatora MPCS jest macierzą utworzoną z dwóch pierwszych wierszy macierzy K.
Macierz K wyznaczono ze wzoru:
K = (MT*M+Λ)−1 * MT , gdzie:
Λ− jest to macierz diagonalna zbudowana ze współczynników kar λ
M − jest macierzą dynamiczną.
Macierz dynamiczną M zbudowano z macierzy odpowiedzi skokowych S powieloną na resztę kolumn macierzy M o wymiarze N x Nu z przesunięciem o wektor 1,1.
Wektor niemierzalnych zakłóceń stanu v(k) wyznaczono z zależności:
x(k) – stan modelu z aktualnej chwili czasowej
x(k-1) – stan modelu z poprzedniej chwili czasowej
u(k-1) sterowanie z poprzedniej chwili czasowej
Błąd modelowania zdefiniowany jest jako różnica między wartością aktualną, a odpowiedzią swobodną modelu i można go wyznaczyć z zależności:
Do regulacji z zastosowaniem regulatora MPCS przyjęto parametry:
horyzontu predykcji N=40
horyzontu sterowania Nu=12
współczynnika kar λ=1
W połowie symulacji po ustabilizowaniu się wartości stężenia oraz temperatury w punkcie pracy
Ca = 260[mol/m3] oraz T=393,9[K] podano na wejście układu skok wartości sterowania stężenia dopływającej substancji Ca=300[mol/m3]. Wartość zadaną temperatury ustawiono na poziomie wartości punktu pracy równej T=393.9[K]. Wartość stężenia Ca oraz temperatury T przedstawiają poniższe charakterystyki.
Stężenie Ca substancji ustawieniu wartości zadanej Ca=300[mol/m3]
Temperatura T substancji ustawieniu wartości zadanej w punkcie pracy T=393.9[K]
Na podstawie otrzymanych charakterystyk stwierdzono, że stężenie substancji stabilizuje się na zadanym poziomie Ca=300[mol/m3] oraz temperatura w punkcie pracy T=393,9[K].
W kolejnym teście w połowie symulacji po ustabilizowaniu się wartości stężenia oraz temperatury w punkcie pracy Ca = 260[mol/m3] oraz T=393,9[K] podano na wejście układu skok wartości zadanej temperatury T=410[K]. Wartość zadaną stężenia ustawiono na poziomie wartości punktu pracy równej Ca=260[mol/m3]. Wartość stężenia Ca oraz temperatury T przedstawiają poniższe charakterystyki.
Stężenie Ca substancji po ustawieniu wartości zadanej w punkcie pracy Ca=260[mol/m3]
Temperatura T substancji po ustawieniu wartości zadanej T=410[K]
Na podstawie otrzymanych charakterystyk stwierdzono, że stężenie stabilizuje się w punkcie pracy Ca=260[ml/m3] oraz temperatura na zadanym poziomieT=410[K].
Zauważono, że stężenie substancji w reaktorze w szybciej osiąga wartość zadaną niż temperatura.
Stwierdzono, że regulator MPCS reaguje znacznie szybciej niż badany regulator PID, który stabilizował charakterystyki stężenia Ca oraz temperatury T dopiero po ok. 10-20tys. iteracji, podczas gdy liczba iteracji niezbędnych do ustabilizowania reakcji przy pomocy regulatora MPCS oscylowała w granicach 1000.
W kolejnym teście ustawiono wartość zadaną stężenia Ca=300[mol/m3] oraz temperatury T=410[K]. Na sterowanie nałożono ograniczenia maksymalnej amplitudy sterowania. Ograniczenie dolne zmiany sterowania stężenia dopływającej substancji Cain ustawiono na Cain = 0, natomiast ograniczenie górne na Cain = 10. Ograniczenie dolne zmiany sterowania przepływem Fc ustawiono na Fc = 0, natomiast ograniczenie górne na Fc = 0, 2.
Stężenie Ca substancji po ustawieniu wartości zadanej Ca=300[mol/m3] oraz T=410[K]
z ograniczeniem amplitudy sygnału sterującego
Temperatura T substancji po ustawieniu wartości zadanej Ca=300[mol/m3] oraz T=410[K]
z ograniczeniem amplitudy sygnału sterującego
Na podstawie przeprowadzonych symulacji stwierdzono, że czas regulacji obiektu sterowanego regulatorem MPCS z ograniczeniami amplitudy sterowania uległ znacznemu wydłużeniu w stosunku do obiektu z regulatorem bez ograniczeń. Zauważono także, że ograniczenia miały większy wpływ na czas stabilizacji temperatury niż na czas stabilizacji stężenia roztworu w reaktorze. Zastosowanie ograniczeń amplitudy sygnału sterującego pozwoliło na zmniejszenie przeregulowań występujących w odpowiedzi skokowej układu.
Kod analitycznego regulatora MPCS:
Matlab |
---|
clear all; clc; close all; kk=3000; %koniec symulacji %% Regulator MPCS gzad(1:kk/2)=260; gzad(kk/2+1:kk)=300; hzad(1:kk/2)=393.9; hzad(kk/2+1:kk)=420; A=[0.1937 -7.7839; 0.0725 1.3079]; B=[0.0601 2.5151; 0.0039 -0.7245]; C=[1 0; 0 1]; N=40; % Horyzont predykcji Nu=12; % Horyzont sterowania lambda = 1; % wspolczynnik kar powinien byś [M, CtAt, CtV] = MPCSmatrices(A,B,C,N,Nu); %funkcja Tatjewskiego liczy M,CtAt,CtV LAMBDA = lambda*eye(2*Nu); % macierz NuxNu K = inv(M'*M+LAMBDA)*M'; % dla zwyklego MPC %% _ _ _ _ _ _ _S Y M U L A C J A y(1,1:kk)=260; y(2,1:kk)=393.9; u(1,1:kk)=2000; u(2,1:kk)=15; x(1,1:kk)=260; x(2,1:kk)=393.9; for i=2:kk % wartość zadana na horyzoncie N for r=1:N if i+r>=kk yzad(2*r-1,1)=gzad(kk); yzad(2*r,1)=hzad(kk); else yzad(2*r-1,1)=gzad(r+i); yzad(2*r,1)=hzad(r+i); end end v(1:2,i)=x(1:2,i)-(A*x(1:2,i-1)+B*u(1:2,i-1)); du(1:2,i) = K(1:2,:)*(yzad - CtAt*x(:,i) - CtV*B*u(:,i-1) - CtV*v(:,i)); u(1:2,i) = u(1:2,i-1)+du(1:2,i); d(1:2,i)=C*x(1:2,i)-C*(A*x(1:2,i-1)+B*u(1:2,i-1)+v(1:2,i-1)); % if u(1,i)>4000 % u(1,i) = 4000; % end % if u(2,i)>30 % u(2,i) = 30; % end % if u(1,i)<1 % u(1,i) = 1; % end % if u(2,i)<1 % u(2,i) = 1; % end % OBIEKT x(1:2,i+1)=A*x(1:2,i)+B*u(1:2,i); y(1:2,i)=C*x(1:2,i); end figure(1) plot(y(1,:),'b'); grid on; title('Ca'); figure(2) plot(y(2,:),'b'); grid on; title('T'); % figure(1) % plot(y(1,kk/2-80:kk),'b'); % grid on; % title('Ca'); % % figure(2) % plot(y(2,kk/2-80:kk),'b'); % grid on; % title('T'); |
Macierz liczona OFF-LINE (macierz statyczna).
Matlab |
---|
[q h] = size(M); lambda = 0.01; % wspolczynnik kar psi = 1; % posejdon PSI = psi*eye(q); % posejdon macierz NxN LAMBDA = lambda*eye(h); % macierz NuxNu H = 2*(M'*PSI*M+LAMBDA); |
Macierz liczona OFF-LINE (macierz statyczna).
Macierz J:
Matlab |
---|
% Macierz statyczna J nu=2; % liczba sygna≥ůw sterujĻcych w uk≥adzie J(1:2*Nu,1:2*Nu)=0; % deklaracja pustej macierzy J k=1; %licznik pomocniczy for j=1:1:2*Nu %zmieŮ kolumnÍ for i=k:2:2*Nu %zmieŮ wiersz J(i,j)=1; if i>=2*Nu-1 k=k+1; end end end % Macierz A Ampcs=[-J;J;-M;M]; |
Wszystkie następne macierze łącznie z ”f” wyliczane są dynamicznie (w każdej iteracji symulacji).
Matlab |
---|
f(1:2*Nu,i)=(-2)*(M)'*( yzad-y0(1:2*N,i) ); |
Matlab |
---|
b((2*Nu+2*N)*nu,1:kk)=0; % U(k-1) do wyliczenia Umin Umax czyli wektorka b for r=1:Nu U(2*r-1,i)=du(1,i-1); U(2*r,i)=du(2,i-1); end y0(:,i)=CtAt*x(:,i) - CtV*B*u(:,i-1) - CtV*v(:,i); … b(:,i)=[-Umin(:)+U(:,i-1); Umax(:)-U(:,i-1); -Ymin(:,1)+y0(:,i); Ymax(:,1)-y0(:,i) ]; |
Gdzie:
Korzystamy z ostatniego wzoru na J(k) opuszczając ostatni człon co daje:
Matlab |
---|
% optymalizacja zadania kwadratowego options = optimoptions(@quadprog, 'Display', 'off'); du(:,i)=quadprog(H,f(24,i),Ampcs(208,24),b(208,i),[],[],[],[],[], options); |
Ostateczny kod:
Matlab |
---|
clear all; clc; close all; kk=80000; %koniec symulacji % wartoúci zadane gzad(1:kk/2)=260; gzad(kk/2+1:kk)=300; hzad(1:kk/2)=393.9; hzad(kk/2+1:kk)=420; %% ___________________________ Regulator MPCS _____________________________ A=[0.1937 -7.7839; 0.0725 1.3079]; B=[0.0601 2.5151; 0.0039 -0.7245]; C=[1 0; 0 1]; N=40; % Horyzont predykcji Nu=12; % Horyzont sterowania lambda = 1; % wspolczynnik kar powinien byś [M, CtAt, CtV] = MPCSmatrices(A,B,C,N,Nu); % funkcja Tatjewskiego liczy M,CtAt,CtV LAMBDA = lambda*eye(2*Nu); % macierz NuxNu K = inv(M'*M+LAMBDA)*M'; % dla zwyklego MPC % Macierz H [q h] = size(M); psi = 1; % posejdon PSI = psi*eye(q); % posejdon macierz NxN LAMBDA = lambda*eye(h); % macierz NuxNu H = 2*(M'*PSI*M+LAMBDA); % Macierz statyczna J nu=2; % liczba sygna≥ůw sterujĻcych w uk≥adzie J(1:2*Nu,1:2*Nu)=0; % deklaracja pustej macierzy J k=1; %licznik pomocniczy for j=1:1:2*Nu %zmieŮ kolumnÍ for i=k:2:2*Nu %zmieŮ wiersz J(i,j)=1; if i>=2*Nu-1 k=k+1; end end end % Macierz A Ampcs=[-J;J;-M;M]; % Ograniczenia sterowania Umin(1,1)=1000; Umin(2,1)=10; Umax(1,1)=3000; Umax(2,1)=20; for r=1:Nu Umin(2*r-1,1)=Umin(1,1); Umin(2*r,1)=Umin(2,1); Umax(2*r-1,1)=Umax(1,1); Umax(2*r,1)=Umax(2,1); end % Ograniczenia wyjúś Ymin(1,1)=260-60; Ymin(2,1)=393.3-90; Ymax(1,1)=260+60; Ymax(2,1)=393.3+90; for r=1:N Ymin(2*r-1,1)=Ymin(1,1); Ymin(2*r,1)=Ymin(2,1); Ymax(2*r-1,1)=Ymax(1,1); Ymax(2*r,1)=Ymax(2,1); end %% __________________________ Wartoúci poczĻtkowe _________________________ y(1,1:kk)=260; y(2,1:kk)=393.9; u(1,1:kk)=2000; u(2,1:kk)=15; x(1,1:kk)=260; x(2,1:kk)=393.9; du(1:2*Nu,1:2)=0; U(1:2*Nu,1)=0; f(1:2*Nu,1:kk)=0; y0(1:2*N,1:kk)=0; b((2*Nu+2*N)*nu,1:kk)=0; %% ___________________________ S Y M U L A C J A __________________________ for i=2:kk % wartoúś zadana (yzad) na horyzoncie N for r=1:N if i+r>=kk yzad(2*r-1,1)=gzad(kk); yzad(2*r,1)=hzad(kk); else yzad(2*r-1,1)=gzad(r+i); yzad(2*r,1)=hzad(r+i); end end % U(k-1) do wyliczenia Umin Umax czyli wektorka b for r=1:Nu U(2*r-1,i)=du(1,i-1); U(2*r,i)=du(2,i-1); end v(1:2,i)=x(1:2,i)-(A*x(1:2,i-1)+B*u(1:2,i-1)); y0(:,i)=CtAt*x(:,i) - CtV*B*u(:,i-1) - CtV*v(:,i); f(1:2*Nu,i)=(-2)*(M)'*( yzad-y0(1:2*N,i) ); b(:,i)=[-Umin(:)+U(:,i-1); Umax(:)-U(:,i-1); -Ymin(:,1)+y0(:,i); Ymax(:,1)-y0(:,i) ]; % optymalizacja zadania kwadratowego options = optimoptions(@quadprog, 'Display', 'off'); du(:,i)=quadprog(H,f(24,i),Ampcs(208,24),b(208,i),[],[],[],[],[], options); % du(:,i)=quadprog(H,f(:,i),[],[],[],[],[],[],[], options); %%%%%%%%%%% KASIA LICZENIE b % b = [-Umin + repmat(U(:,k-1)', 1, Nu)'; % Umax - repmat(U(:,k-1)', 1, Nu)'; % -Ymin + CtAt*X(:,k) + CtV*B_d*U(:,k-1) + CtV*v(:,k); %-Ymin + Y0 % Ymax - CtAt*X(:,k) - CtV*B_d*U(:,k-1) - CtV*v(:,k)]; % Ymax - Y0 % % ,gdzie % U(2, 1:2) = Fh0; % U(1, 1:2) = Fc0; % B_d to nasza macierz B % % Umin = (repmat([0 0], 1, Nu))'; % Umax = (repmat([Fc_max Fh_in_max], 1, Nu))'; % dUmax = (repmat([dFc_max dFh_in_max], 1, Nu))'; % Ymin = (repmat([h_min Tout_min], 1,N))'; % Ymax = (repmat([h_max Tout_max], 1,N))'; %%%%%%%%%%% End Kasia u(1:2,i) = u(1:2,i-1)+du(1:2,i); % d(1:2,i)=C*x(1:2,i)-C*(A*x(1:2,i-1)+B*u(1:2,i-1)+v(1:2,i-1)); % OBIEKT x(1:2,i+1)=A*x(1:2,i)+B*u(1:2,i); y(1:2,i)=C*x(1:2,i); end figure(1) plot(y(1,:),'b'); grid on; title('Ca'); figure(2) plot(y(2,:),'b'); grid on; title('T'); |