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

Zadanie

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

Dwupętlowy układ regulacji z regulatorami PID
1. Układ regulacji bez odprzęgania

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/m³] na 300 [mol/m³].

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.

Układ regulacji z odsprzęganiem

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.

Układ regulacji z regulatorem MPCS
1. Regulator analityczny MPCS

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 = (M^T*M+Λ)⁻¹ * M^T , 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 = lambdaeye(2Nu); % 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(2r-1,1)=gzad(kk); yzad(2r,1)=hzad(kk); else yzad(2r-1,1)=gzad(r+i); yzad(2r,1)=hzad(r+i); end end v(1:2,i)=x(1:2,i)-(Ax(1:2,i-1)+Bu(1:2,i-1)); du(1:2,i) = K(1:2,:)(yzad - CtAtx(:,i) - CtVBu(:,i-1) - CtVv(:,i)); u(1:2,i) = u(1:2,i-1)+du(1:2,i); d(1:2,i)=Cx(1:2,i)-C(Ax(1:2,i-1)+Bu(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)=Ax(1:2,i)+Bu(1:2,i); y(1:2,i)=Cx(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');

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

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');

Regulator numeryczny MPCS

Macierz H:

Macierz liczona OFF-LINE (macierz statyczna).

Matlab
[q h] = size(M); lambda = 0.01; % wspolczynnik kar psi = 1; % posejdon PSI = psieye(q); % posejdon macierz NxN LAMBDA = lambdaeye(h); % macierz NuxNu H = 2(M'PSI*M+LAMBDA);

[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 A:

Macierz liczona OFF-LINE (macierz statyczna).

Macierz J:

Matlab
% Macierz statyczna J nu=2; % liczba sygna≥ůw sterujĻcych w uk≥adzie J(1:2Nu,1:2Nu)=0; % deklaracja pustej macierzy J k=1; %licznik pomocniczy for j=1:1:2Nu %zmieŮ kolumnÍ for i=k:2:2Nu %zmieŮ wiersz J(i,j)=1; if i>=2*Nu-1 k=k+1; end end end % Macierz A Ampcs=[-J;J;-M;M];

% 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

% Macierz A

Ampcs=[-J;J;-M;M];

Macierz f:

Wszystkie następne macierze łącznie z ”f” wyliczane są dynamicznie (w każdej iteracji symulacji).

Matlab
f(1:2Nu,i)=(-2)(M)'( yzad-y0(1:2N,i) );

Macierz b:

Matlab
b((2Nu+2N)nu,1:kk)=0; % U(k-1) do wyliczenia Umin Umax czyli wektorka b for r=1:Nu U(2r-1,i)=du(1,i-1); U(2r,i)=du(2,i-1); end y0(:,i)=CtAtx(:,i) - CtVBu(:,i-1) - CtV*v(:,i); … b(:,i)=[-Umin(:)+U(:,i-1); Umax(:)-U(:,i-1); -Ymin(:,1)+y0(:,i); Ymax(:,1)-y0(:,i) ];

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) ];

Macierz J:

Gdzie:

Korzystamy z ostatniego wzoru na J(k) opuszczając ostatni człon co daje:

Zadanie optymalizacji dla zadania kwadratowego:

Matlab
% optymalizacja zadania kwadratowego options = optimoptions(@quadprog, 'Display', 'off'); du(:,i)=quadprog(H,f(24,i),Ampcs(208,24),b(208,i),[],[],[],[],[], options);

% 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 = lambdaeye(2Nu); % macierz NuxNu K = inv(M'M+LAMBDA)M'; % dla zwyklego MPC % Macierz H [q h] = size(M); psi = 1; % posejdon PSI = psieye(q); % posejdon macierz NxN LAMBDA = lambdaeye(h); % macierz NuxNu H = 2(M'PSIM+LAMBDA); % Macierz statyczna J nu=2; % liczba sygna≥ůw sterujĻcych w uk≥adzie J(1:2Nu,1:2Nu)=0; % deklaracja pustej macierzy J k=1; %licznik pomocniczy for j=1:1:2Nu %zmieŮ kolumnÍ for i=k:2:2Nu %zmieŮ wiersz J(i,j)=1; if i>=2Nu-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(2r-1,1)=Umin(1,1); Umin(2r,1)=Umin(2,1); Umax(2r-1,1)=Umax(1,1); Umax(2r,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(2r-1,1)=Ymin(1,1); Ymin(2r,1)=Ymin(2,1); Ymax(2r-1,1)=Ymax(1,1); Ymax(2r,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:2Nu,1:2)=0; U(1:2Nu,1)=0; f(1:2Nu,1:kk)=0; y0(1:2N,1:kk)=0; b((2Nu+2N)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(2r-1,1)=gzad(kk); yzad(2r,1)=hzad(kk); else yzad(2r-1,1)=gzad(r+i); yzad(2r,1)=hzad(r+i); end end % U(k-1) do wyliczenia Umin Umax czyli wektorka b for r=1:Nu U(2r-1,i)=du(1,i-1); U(2r,i)=du(2,i-1); end v(1:2,i)=x(1:2,i)-(Ax(1:2,i-1)+Bu(1:2,i-1)); y0(:,i)=CtAtx(:,i) - CtVBu(:,i-1) - CtVv(:,i); f(1:2Nu,i)=(-2)(M)'( yzad-y0(1:2N,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 + CtAtX(:,k) + CtVB_dU(:,k-1) + CtVv(:,k); %-Ymin + Y0 % Ymax - CtAtX(:,k) - CtVB_dU(:,k-1) - CtVv(:,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)=Cx(1:2,i)-C(Ax(1:2,i-1)+Bu(1:2,i-1)+v(1:2,i-1)); % OBIEKT x(1:2,i+1)=Ax(1:2,i)+Bu(1:2,i); y(1:2,i)=Cx(1:2,i); end figure(1) plot(y(1,:),'b'); grid on; title('Ca'); figure(2) plot(y(2,:),'b'); grid on; title('T');

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

% 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

% 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');

Wnioski

Wyszukiwarka

TAP Projekt Etap 2 Szyszka Terwiński Krężel

Zadanie

Dwupętlowy układ regulacji z regulatorami PID

Układ regulacji bez odprzęgania

Układ regulacji z odsprzęganiem

Układ regulacji z regulatorem MPCS

Regulator analityczny MPCS

Regulator numeryczny MPCS

Macierz H:

Macierz A:

Macierz f:

Macierz b:

Macierz J:

Zadanie optymalizacji dla zadania kwadratowego:

Wnioski