Jeżeli do predykcji stosuje się model liniowy, wyjściowe sygnały prognozowane \(\hat{y}(k+1|k),\ldots,\hat{y}(k+N|k)\) są liniowymi funkcjami przyszłych sygnałów sterujących - zmiennych decyzyjnych algorytmu regulacji predykcyjnej. W konsekwencji, problem optymalizacji algorytmu jest zadaniem optymalizacji kwadratowej (kwadratowa funkcja celu, liniowe ograniczenia). Minimum globalne takiego zadania może być łatwo wyznaczone przy użyciu algorytmu ograniczeń aktywnych lub algorytmu punktu wewnętrznego w skończonej liczbie iteracji, można przewidzieć czas obliczeń. Algorytmy regulacji predykcyjnej bazujące na modelach liniowych (zwane krótko algorytmami liniowymi), przede wszystkim Dynamic Matrix Control (DMC) [CR79] oraz Generalized Predictive Control (GPC) [CMT87] odniosły niekwestionowany sukces w praktyce, wywierając tym samym istotny wpływ na kierunki prac badawczych [T07, T02]. W praktyce pracują one znacznie lepiej niż stosowane uprzednio algorytmy typu PID (najczęściej jednopętlowe).
Właściwości bardzo wielu procesów technologicznych są zwykle nieliniowe. Zastosowanie algorytmu liniowego do nieliniowego procesu jest często możliwe, ale otrzymana jakość regulacji nie jest wówczas zadowalająca, szczególnie przy następujących szybko i w szerokim zakresie zmianach punktu pracy. Jest to spowodowane tym, że model liniowy jedynie w ograniczonym zakresie opisuje zachowanie procesu nieliniowego. W przypadku procesów nieliniowych naturalnym rozwiązaniem jest zastosowanie nieliniowych algorytmów regulacji predykcyjnej, w których do predykcji stosuje się modele nieliniowe. Dwa zasadnicze problemy, przed którymi stoi projektant algorytmu nieliniowego, są następujące:
Niech nieliniowy model dyskretny procesu dynamicznego ma ogólną postać \begin{equation} y(k)=f(\boldsymbol{x}(k))=f(u(k-\tau),\ldots,u(k-n_{\mathrm{B}}),y(k-1),\ldots,y(k-n_{\mathrm{A}})) \label{w_modnl_ogolnie} \end{equation} gdzie \(\tau\) jest opóźnieniem, liczby naturalne \(n_{\mathrm{A}}\) i \(n_{\mathrm{B}}\) określają rząd dynamiki modelu, natomiast \(f\colon \mathbb{R}^{n_{\mathrm{A}}+n_{\mathrm{B}}-\tau+1}\rightarrow\mathbb{R}\) jest nieliniową funkcją charakteryzującą model. Stosując model rekurencyjnie, kolejne wyjściowe sygnały prognozowane na horyzoncie predykcji (\(p=1,\ldots,N\)) można wyrazić następująco \begin{align} \hat{y}(k+p|k) =f(&\underbrace{u(k-\tau+p|k),\ldots,u(k|k)}_{I_{\mathrm{uf}}(p)},% \underbrace{u(k-1),\ldots,u(k-n_{\mathrm{B}}+p)}_{I_\mathrm{u}-I_{\mathrm{uf}}(p)},\nonumber\\ & \underbrace{\hat{y}(k-1+p|k),\ldots,\hat{y}(k+1|k)}_{I_{\mathrm{yf}}(p)}, \underbrace{y(k),\ldots,y(k-n_{\mathrm{A}}+p)}_{n_{\mathrm{A}}-I_{\mathrm{yf}}(p)})\nonumber\\ &+d(k)\label{w_ypredkpk_f} \end{align} Sygnały prognozowane są funkcją: \(I_{\mathrm{uf}}(p)=\max(\min(p-\tau+1,I_{\mathrm{u}}),0)\) (gdzie \(I_{\mathrm{u}}=n_{\mathrm{B}}-\tau+1\)) przyszłych wartości sygnału sterującego (czyli zmiennych decyzyjnych algorytmu), \(I_{\mathrm{u}}-I_{\mathrm{uf}}(p)\) wartości sygnału sterującego, które zostały wykorzystane do sterowania w poprzednich iteracjach, \(I_{\mathrm{yf}}(p)=\min(p-1,n_{\mathrm{A}})\) przewidywanych wartości sygnału wyjściowego oraz \(n_{\mathrm{A}}-I_{\mathrm{yf}}(p)\) wartości sygnału wyjściowego procesu zmierzonych w poprzednich iteracjach, \(d(k)\) jest oszacowaniem niemierzalnego zakłócenia działąjącego na wyjście procesu. Warto podkreślić, że w wyniku zastosowania do predykcji modelu nieliniowego, sygnały prognozowane są nieliniowymi funkcjami aktualnie obliczanymi przyszłymi sygnałami sterującymi. W rezultacie otrzymuje się nieliniowe zadanie optymalizacji, które musi być rozwiązywane cyklicznie, na bieżąco (ang. on-line) w każdej iteracji algorytmu. Nieliniowa optymalizacja jest nie tylko złożona obliczeniowo, ale też może okazać się zawodna, gdyż pojawia się problem minimów lokalnych. Oczywiście, czas obliczeń nie może być zagwarantowany.
Praktycznym rozwiązaniem jest cykliczna, bieżąca linearyzacja modelu nieliniowego dla aktualnego punktu pracy (ang. successive linearisation). Do predykcji stosuje się model liniowy o parametrach zmiennych w czasie, co prowadzi do kwadratowego zadania optymalizacji. Algorytmy takie są oczywiście suboptymalne, gdyż model zlinearyzowany jest jedynie pewnym przybliżeniem modelu nieliniowego, a tym samym procesu. Suboptymalne algorytmy regulacji predykcyjnej są przez praktyków cenione i stosowane. M. Morari i J. Lee stwierdzają: „Linearyzacja jest jedyną metodą, która znalazła szersze zastosowanie w przemyśle, pomijając projekty demonstracyjne. Dla przemysłu potrzeba rozwiązywania nieliniowych zadań dynamicznych w czasie rzeczywistym musi być jasno uzasadniona, a jak dotąd nie ma przykładów wskazujących taką potrzebę w sposób przekonujący” [ML99]. Przedstawione poniżej algorytmy suboptymalne są uniwersalne, można w nich zastosować różne klasy modeli. Jedynym warunkiem jest różniczkowalność modelu.
W najprostszym algorytmie suboptymalnym nieliniowy model procesu \eqref{w_modnl_ogolnie} jest cyklicznie linearyzowany w aktualnym punkcie pracy [T07, T02]. Korzystając z rozwinięcia funkcji w szereg Taylora, model zlinearyzowany, czyli aproksymacja liniowa nieliniowej funkcji \(f(\boldsymbol{x}(k))\colon \mathbb{R}^{n_{\mathrm{A}}+n_{\mathrm{B}}-\tau+1}\rightarrow\mathbb{R}\) ma postać \begin{align} y(k)&=f(\bar{\boldsymbol{x}}(k))+\sum_{l=\tau}^{n_{\mathrm{B}}}b_l(\bar{\boldsymbol{x}}(k))(u(k-l)-\bar{u}(k-l))\nonumber\\ &\quad -\sum_{l=1}^{n_{\mathrm{A}}}a_l(\bar{\boldsymbol{x}}(k))(y(k-l)-\bar{y}(k-l)) %\nonumber \label{w_mod_zlin} \end{align} przy czym wektor pomiarów \(\bar{\boldsymbol{x}}(k)\) definiuje aktualny punkt pracy, natomiast współczynniki modelu zlinearyzowanego obliczane są jako pochodne modelu nieliniowego względem jego argumentów w aktualnym punkcie pracy \begin{equation} a_l(k)=-{\left. \frac{\partial f(\boldsymbol{x}(k))}{\partial y(k-l)}\right| _{\boldsymbol{x}(k)=\bar{\boldsymbol{x}}(k)}}, \ b_l(k)={\left. \frac{\partial f(\boldsymbol{x}(k))}{\partial u(k-l)}\right| _{\boldsymbol{x}(k)=\bar{\boldsymbol{x}}(k)}}\nonumber \end{equation} Dzięki linearyzacji otrzymuje się liniowe równanie predykcji \begin{equation} \hat{\boldsymbol{y}}(k)= \left[ \begin{array} {c}% \hat{y}(k+1)\\ \vdots\\ \hat{y}(k+N) \end{array} \right]=\boldsymbol{G}(k)\triangle\boldsymbol{u}(k) + \boldsymbol{y}^{0}(k) %\nonumber \label{w_y_Gdu_y0} \end{equation} gdzie macierz zależnych od punktu pracy współczynników odpowiedzi skokowej modelu zlinearyzowanego ma postać \begin{equation} \boldsymbol{G}(k)=\left[ \begin{array} [c]{cccc}% s_{1}(k) & 0 & \ldots & 0\\[0.25mm] s_{2}(k) & s_{1}(k) & \ldots & 0\\[0.25mm] \vdots & \vdots & \ddots & \vdots\\ s_{N}(k) & s_{N-1}(k) & \ldots & s_{N-N_{\mathrm{u}}+1}(k)% \end{array} \right] \nonumber \end{equation} natomiast \(\boldsymbol{y}^{0}(k)\) jest wektorem trajektorii swobodnej, która opisuje reakcję procesu na przeszłość. Trajektoria ta jest liczona przy wykorzystaniu modelu zlinearyzowanego \eqref{w_mod_zlin}.
Stosując suboptymalne równanie predykcji \eqref{w_y_Gdu_y0}, zamiast nieliniowego zadania optymalizacji, w algorytmie SL otrzymuje się zadanie optymalizacji kwadratowej \begin{align} &\begin{array}{l} \min\limits_{\triangle \boldsymbol{u}(k)}\Big\{J(k)=\left\| \boldsymbol{y} ^{\mathrm{zad}}(k)-\boldsymbol{G}(k)\triangle\boldsymbol{u}(k)-\boldsymbol{y}^{0}(k)\right\|^{2}_{\boldsymbol{M}} + \left\| \triangle\boldsymbol{u}(k)\right\| ^{2}_{\boldsymbol{\Lambda}} \Big\} \end{array}\nonumber\\ &\text{przy ograniczeniach}\label{w_zad_opt_sl}\\ &\boldsymbol{u}^{\min}\leq\boldsymbol{J}\triangle\boldsymbol{u}(k)+\boldsymbol{u}(k-1) \leq\boldsymbol{u}^{\max}\nonumber\\ &-\triangle\boldsymbol{u}^{\max}\leq\triangle\boldsymbol{u}(k)\leq\triangle\boldsymbol{u}^{\max}\nonumber\\ &\boldsymbol{y}^{\min}\leq\boldsymbol{G}(k)\triangle\boldsymbol{u}(k)+\boldsymbol{y}^{0}(k)\leq\boldsymbol{y}^{\max}\nonumber \end{align} gdzie \(\boldsymbol{M}=\mathrm{diag}(\boldsymbol{M}_1,\ldots,\boldsymbol{M}_N)\), \(\boldsymbol{\Lambda}=\mathrm{diag}(\boldsymbol{\Lambda}_0,\ldots,\boldsymbol{\Lambda}_{N_{\mathrm{u}}-1})\), \(\boldsymbol{J}\) jest macierzą dolną trójkątną zawierającą elementy jednostkowe.
Struktura algorytmu SL została przedstawiona na rys. 2. W każdej iteracji algorytmu nieliniowy model procesu \eqref{w_modnl_ogolnie} jest cyklicznie linearyzowany, na jego pdstawie obliczana jest macierz współczynników odpowiedzi skokowej \(\boldsymbol{G}(k)\) oraz trajektoria swobodna \(\boldsymbol{y}^{0}(k)\), a następnie rozwiązywane jest zadanie optymalizacji kwadratowej \eqref{w_zad_opt_sl} i pierwszy element wektora zmiennych decyzyjnych jest stosowany do sterowania. W kolejnej iteracji, po aktualizacji pomiarów, horyzont zostaje przesunięty jeden krok do przodu i cała procedura obliczeniowa zostaje powtórzona.
Istnieje prosta metoda zwiększenia jakości regulacji algorytmu SL. Można mianowicie założyć, że trajektoria swobodna \(\boldsymbol{y}^{0}(k)\) obliczana jest nie przy użyciu modelu zlinearyzowanego w aktualnym punkcie pracy \eqref{w_mod_zlin}, lecz oryginalnego modelu nieliniowego \eqref{w_modnl_ogolnie} [T07, T02]. Struktura algorytmu NPL została przedstawiona na rys. 3. Warto zauważyć, że ogólna sieć działań algorytmów SL i NPL jest taka sama, rozwiązuje się to samo zadanie optymalizacji \eqref{w_zad_opt_sl}, jedyna różnica polega na użyciu nieliniowej trajektorii swobodnej. Z punktu widzenia złożoności obliczeniowej jest to praktycznie bez znaczenia, ponieważ decydujący wpływ na złożoność obliczeniową algorytmu ma procedura optymalizacji kwadratowej, natomiast, jak wykazują przeprowadzone badania [T07, T02], zastosowanie nieliniowej trajektorii swobodnej znacznie poprawia jakość regulacji w porównaniu z algorytmem SL.
Cechą szczególną suboptymalnych algorytmów SL i NPL jest bieżąca, cykliczna linearyzacja modelu nieliniowego w aktualnym punkcie pracy i wykorzystanie tego samego modelu zlinearyzowanego do predykcji na całym horyzoncie predykcji. Jeżeli linearyzowany (nieliniowy) model jest dokładny, natomiast algorytmy suboptymalne SL i NPL pracują znacznie gorzej niż algorytm z nieliniową optymalizacją powtarzaną w każdej iteracji, oczywiste jest, że wykorzystująca zlinearyzowany model predykcja suboptymalna \eqref{w_y_Gdu_y0} jest znacząco różna od dokładnej predykcji nieliniowej \eqref{w_ypredkpk_f}. Jest to spowodowane tym, że przy dużych zmianach punktu pracy lub też gdy występują silne zakłócenia, lokalne przybliżenie liniowe może aproksymować właściwości modelu nieliniowego (a tym samym rzeczywistego procesu) jedynie na początku horyzontu predykcji. Przy długich horyzontach różnice między trajektorią suboptymalną a trajektorią nieliniową (i trajektorią procesu) mogą być znaczące.
Uwzględniając opisane powyżej niedoskonałości linearyzacji w punkcie pracy, naturalne jest zastosowanie linearyzacji wokół pewnej przyszłej trajektorii sygnałów wejściowych [Ł14a, Ł13] \begin{equation} \boldsymbol{u}^{\mathrm{traj}}(k)=\left[ \begin{array}{c} u^{\mathrm{traj}}(k|k)\\ \vdots\\ u^{\mathrm{traj}}(k+N_{\mathrm{u}}-1|k) \end{array} \right] \nonumber \end{equation} Założonej trajektorii wejściowej \(\boldsymbol{u}^{\mathrm{traj}}(k)\) odpowiada prognozowana trajektoria wyjściowa (może być ona obliczona przy wykorzystaniu nieliniowego modelu \eqref{w_modnl_ogolnie}) \begin{equation} \hat{\boldsymbol{y}}^{\mathrm{traj}}(k)=\left[ \begin{array}{c} y^{\mathrm{traj}}(k+1|k)\\ \vdots\\ y^{\mathrm{traj}}(k+N|k) \end{array} \right] \nonumber \end{equation} Korzystając z rozwinięcia funkcji w szereg Taylora, aproksymacja liniowa nieliniowej trajektorii prognozowanej sygnałów wyjściowych \(\hat{\boldsymbol{y}}(k)\), czyli funkcji \(\hat{\boldsymbol{y}}(\boldsymbol{u}(k))\colon \mathbb{R}^{N_{\mathrm{u}}}\rightarrow\mathbb{R}^{N}\), ma postać \begin{equation} \hat{\boldsymbol{y}}(k)=\hat{\boldsymbol{y}}^{\mathrm{traj}}(k)+\boldsymbol{H}(k)(\boldsymbol{u}(k)-\boldsymbol{u}^{\mathrm{traj}}(k))\nonumber \end{equation} gdzie macierz pochodnych cząstkowych prognozowanej trajektorii wyjściowej względem przyszłej trajektorii sterującej jest obliczana dla założonej trajektorii \(\boldsymbol{u}^{\mathrm{traj}}(k)\) \begin{equation} \boldsymbol{H}(k)={\left. \dfrac{\mathrm{d}\hat{\boldsymbol{y}}(k)}{\mathrm{d}\boldsymbol{u}(k)}\right| _{{\begin{smallmatrix} \hat{\boldsymbol{y}}(k)=\hat{\boldsymbol{y}}^{\mathrm{traj}}(k)\\ \boldsymbol{u}(k)=\boldsymbol{u}^{\mathrm{traj}}(k) \end{smallmatrix}}}} =\left[ \begin{array} [c]{ccc}% \dfrac{\partial \hat{y}^{\mathrm{traj}}(k+1|k)}{\partial u^{\mathrm{traj}}(k|k)} & \cdots & \dfrac{\partial \hat{y}^{\mathrm{traj}}(k+1|k)}{\partial u^{\mathrm{traj}}(k+N_{\mathrm{u}}-1|k)}\\ \vdots & \ddots & \vdots\\ \dfrac{\partial \hat{y}^{\mathrm{traj}}(k+N|k)}{\partial u^{\mathrm{traj}}(k|k)} & \cdots & \dfrac{\partial \hat{y}^{\mathrm{traj}}(k+N|k)}{\partial u^{\mathrm{traj}}(k+N_{\mathrm{u}}-1|k)}% \end{array} \right] \nonumber \end{equation} Liniową aproksymację nieliniowej trajektorii prognozowanej sygnałów wyjściowych można przedstawić jako funkcję przyszłych przyrostów sygnałów sterujących \(\triangle \boldsymbol{u}(k)\), czyli zmiennych decyzyjnych algorytmu regulacji \begin{equation} \hat{\boldsymbol{y}}(k)=\boldsymbol{H}(k)\boldsymbol{J}\triangle \boldsymbol{u}(k)+\hat{\boldsymbol{y}}^{\mathrm{traj}}(k)+\boldsymbol{H}(k)(\boldsymbol{u}(k-1)-\boldsymbol{u}^{\mathrm{traj}}(k))\label{w_zasada_superpozycji_nplt} \end{equation} Dzięki linearyzacji wyjściowa trajektoria prognozowana jest liniową funkcją zmiennych decyzyjnych algorytmu, analogicznie jak w przypadku algorytmów SL i NPL (równanie \eqref{w_y_Gdu_y0}).
Stosując suboptymalne równanie predykcji \eqref{w_zasada_superpozycji_nplt}, zamiast nieliniowego zadania optymalizacji, w algorytmie NPLT otrzymuje się zadanie optymalizacji kwadratowej \begin{align} &\begin{array} {l} \min\limits_{\triangle \boldsymbol{u}(k) }\Big\{J(k)=\big\| \boldsymbol{y} ^{\mathrm{zad}}(k)-\boldsymbol{H}(k)\boldsymbol{J}\triangle \boldsymbol{u}(k)-\hat{\boldsymbol{y}}^{\mathrm{traj}}(k)\\ \ \ \, \quad \quad \quad \quad \qquad - \: \boldsymbol{H}(k)(\boldsymbol{u}(k-1)-\boldsymbol{u}^{\mathrm{traj}}(k))\big\|^{2}_{\boldsymbol{M}} + \left\|\triangle\boldsymbol{u}(k)\right\| ^{2}_{\boldsymbol{\Lambda}} \Big\} \end{array}\nonumber\\ &\text{przy ograniczeniach}\label{w_zad_opt_nplt}\\ &\boldsymbol{u}^{\min}\leq\boldsymbol{J}\triangle\boldsymbol{u}(k)+\boldsymbol{u}(k-1) \leq\boldsymbol{u}^{\max}\nonumber\\ &-\triangle\boldsymbol{u}^{\max}\leq\triangle\boldsymbol{u}(k)\leq\triangle\boldsymbol{u}^{\max}\nonumber\\ &\boldsymbol{y}^{\min}\leq \boldsymbol{H}(k)\boldsymbol{J}\triangle \boldsymbol{u}(k)+\hat{\boldsymbol{y}}^{\mathrm{traj}}(k)+\boldsymbol{H}(k)(\boldsymbol{u}(k-1)-\boldsymbol{u}^{\mathrm{traj}}(k)) \leq\boldsymbol{y}^{\max}\nonumber \end{align}
Struktura algorytmu NPLT została przedstawiona na rys. 4. W każdej iteracji algorytmu, na podstawie nieliniowego modelu procesu, wyznacza się prognozowaną trajektorię wyjściową \(\hat{\boldsymbol{y}}^{\mathrm{traj}}(k)\) odpowiadającą założonej przyszłej trajektorii sterującej \(\boldsymbol{u}^{\mathrm{traj}}(k)\) oraz wyznacza się aproksymację liniową trajektorii prognozowanej \(\hat{\boldsymbol{y}}(k)\) wokół trajektorii \(\boldsymbol{u}^{\mathrm{traj}}(k)\) daną wzorem \eqref{w_zasada_superpozycji_nplt}, tzn. macierz pochodnych cząstkowych trajektorii prognozowanej \(\boldsymbol{H}(k)\). Następnie rozwiązywane jest zadanie optymalizacji kwadratowej \eqref{w_zad_opt_nplt} i pierwszy element wektora zmiennych decyzyjnych jest stosowany do sterowania. W kolejnej iteracji, po aktualizacji pomiarów, horyzont zostaje przesunięty jeden krok do przodu i cała procedura obliczeniowa zostaje powtórzona.
Sprawą zasadniczą jest wybór przyszłej trajektorii sygnałów wejściowych \(\boldsymbol{u}^{\mathrm{traj}}(k)\), wokół której wykonywana jest linearyzacja prognozowanej trajektorii wyjściowej \(\hat{\boldsymbol{y}}(k)\). Ideałem jest linearyzacja dla trajektorii optymalnej, obliczanej w aktualnej iteracji algorytmu, ale nie jest ona oczywiście znana przed rozwiązaniem zadania optymalizacji. Najprostszym rozwiązaniem jest wykorzystanie do linearyzacji trajektorii zdefiniowanej przez sygnał sterujący wyznaczony (i zastosowany do sterowania) w poprzedniej iteracji \begin{equation} \boldsymbol{u}^{\mathrm{traj}}(k)=\boldsymbol{u}(k-1)=\left[ \begin{array} {c}% u(k-1)\\ \vdots\\ u(k-1) \end{array} \right] \nonumber \end{equation} Trajektoria wyjściowa \(\hat{\boldsymbol{y}}^{\mathrm{traj}}(k)\) jest wówczas nieliniową trajektorią swobodną. Algorytm NPLT z linearyzacją wokół trajektorii \(\boldsymbol{u}(k-1)\) jest oznaczany jako \(\mathrm{NPLT}_{\boldsymbol{u}(k-1)}\). Można również zastosować linearyzację wokół \(N_{\mathrm{u}}-1\) ostatnich elementów trajektorii optymalnej wyznaczonej w poprzedniej iteracji algorytmu \begin{equation} \boldsymbol{u}^{\mathrm{traj}}(k)=\boldsymbol{u}(k|k-1)=\left[ \begin{array} {c}% u(k|k-1)\\ \vdots\\ u(k+N_{\mathrm{u}}-3|k-1)\\[0.5mm] u(k+N_{\mathrm{u}}-2|k-1)\\[0.5mm] u(k+N_{\mathrm{u}}-2|k-1) \end{array} \right] \nonumber \end{equation} Druga wersja algorytmu NPLT oznaczona jest jako \(\mathrm{NPLT}_{\boldsymbol{u}(k|k-1)}\).
W algorytmie NPLT przyjmuje się arbitralnie pewną trajektorię zmian przyszłych sygnałów sterujących na horyzoncie sterowania \(\boldsymbol{u}^{\mathrm{traj}}(k)\) i dokonuje się linearyzacji wyjściowej trajektorii prognozowanej \(\hat{y}(k)\) wokół tej trajektorii. Algorytm NPLT jest polecany wówczas, gdy jakość regulacji klasycznego algorytmu NPL z linearyzacją modelu w aktualnym punkcie pracy jest niezadowalająca. W każdej iteracji \(k\) algorytmu linearyzacja wokół przyszłej trajektorii sterującej oraz obliczenie zmiennych decyzyjnych mogą być wykonywane kilkukrotnie (iteracyjnie) [Ł14a, Ł13, Ł11c]. Intuicyjnie, w porównaniu z algorytmem NPLT, powinno to doprowadzić do zwiększenia jakości predykcji, a tym samym poprawić jakość regulacji.
Niech \(t\) będzie indeksem iteracji wewnętrznych (\(t=1,2,3,\ldots\)). Korzystając z rozwinięcia funkcji w szereg Taylora, aproksymacja liniowa nieliniowej trajektorii prognozowanej sygnałów wyjściowych \(\hat{\boldsymbol{y}}^t(k)\), czyli funkcji \(\hat{\boldsymbol{y}}^t(\boldsymbol{u}^t(k))\colon \mathbb{R}^{N_{\mathrm{u}}}\rightarrow\mathbb{R}^{N}\), ma postać \begin{equation} \hat{\boldsymbol{y}}^t(k)=\hat{\boldsymbol{y}}^{t-1}(k)+\boldsymbol{H}^t(k)(\boldsymbol{u}^t(k)-\boldsymbol{u}^{t-1}(k))\label{w_zasada_superpozycji_nplptp0} \end{equation} gdzie macierz pochodnych cząstkowych prognozowanej trajektorii wyjściowej względem przyszłej trajektorii sterującej jest obliczana dla trajektorii wejściowej \(\boldsymbol{u}^{t-1}(k)\) wyznaczonej w poprzedniej iteracji wewnętrznej \begin{equation} \boldsymbol{H}^t(k)={\left. \dfrac{\mathrm{d}\hat{\boldsymbol{y}}(k)}{\mathrm{d}\boldsymbol{u}(k)}\right| _{{\begin{smallmatrix} \hat{\boldsymbol{y}}(k)=\hat{\boldsymbol{y}}^{t-1}(k)\\ \boldsymbol{u}(k)=\boldsymbol{u}^{t-1}(k) \end{smallmatrix}}}} =\left[ \begin{array} [c]{ccc}% \dfrac{\partial \hat{y}^{t-1}(k+1|k)}{\partial u^{t-1}(k|k)} & \cdots & \dfrac{\partial \hat{y}^{t-1}(k+1|k)}{\partial u^{t-1}(k+N_{\mathrm{u}}-1|k)}\\ \vdots & \ddots & \vdots\\ \dfrac{\partial \hat{y}^{t-1}(k+N|k)}{\partial u^{t-1}(k|k)} & \cdots & \dfrac{\partial \hat{y}^{t-1}(k+N|k)}{\partial u^{t-1}(k+N_{\mathrm{u}}-1|k)} \end{array} \right]\nonumber \end{equation}
Stosując suboptymalne równanie predykcji \eqref{w_zasada_superpozycji_nplptp0}, zamiast nieliniowego zadania optymalizacji, w algorytmie NPLTP otrzymuje się zadanie optymalizacji kwadratowej \begin{align} &\begin{array}{l} \min\limits_{\triangle \boldsymbol{u}^t(k)}\Big\{J(k)=\big\| \boldsymbol{y} ^{\mathrm{zad}}(k)-\boldsymbol{H}^t(k)\boldsymbol{J}\triangle \boldsymbol{u}^t(k)-\hat{\boldsymbol{y}}^{t-1}(k)\\ \ \quad \quad \qquad \qquad \ \, - \: \boldsymbol{H}^t(k)(\boldsymbol{u}(k-1)-\boldsymbol{u}^{t-1}(k))\big\|^{2}_{\boldsymbol{M}} + \left\|\triangle\boldsymbol{u}^t(k)\right\| ^{2}_{\boldsymbol{\Lambda}} \Big\} \end{array}\nonumber\\ &\text{przy ograniczeniach}\label{w_zad_opt_npltp}\\ &\boldsymbol{u}^{\min}\leq\boldsymbol{J}\triangle\boldsymbol{u}^t(k)+\boldsymbol{u}(k-1) \leq\boldsymbol{u}^{\max}\nonumber\\ &-\triangle\boldsymbol{u}^{\max}\leq\triangle\boldsymbol{u}^t(k)\leq\triangle\boldsymbol{u}^{\max}\nonumber\\ &\boldsymbol{y}^{\min}\leq \boldsymbol{H}^t(k)\boldsymbol{J}\triangle \boldsymbol{u}^t(k)+\hat{\boldsymbol{y}}^{t-1}(k) +\boldsymbol{H}^t(k)(\boldsymbol{u}(k-1)-\boldsymbol{u}^{t-1}(k)) \leq\boldsymbol{y}^{\max}\nonumber \end{align}
Struktura algorytmu NPLTP została przedstawiona na rys. 5. W każdej iteracji wewnętrznej \(t\) algorytmu, na podstawie nieliniowego modelu procesu, wyznacza się prognozowaną trajektorię wyjściową \(\hat{\boldsymbol{y}}^{t-1}(k)\) odpowiadającą przyszłej trajektorii sterującej \(\boldsymbol{u}^{t-1}(k)\) obliczonej w poprzedniej iteracji wewnętrznej oraz wyznacza się aproksymację liniową trajektorii prognozowanej \(\hat{\boldsymbol{y}}^t(k)\) wokół trajektorii trajektorii wejściowej \(\boldsymbol{u}^{t-1}(k)\) daną wzorem \eqref{w_zasada_superpozycji_nplptp0}, tzn. macierz pochodnych cząstkowych trajektorii prognozowanej \(\boldsymbol{H}^t(k)\). Następnie rozwiązywane jest zadanie optymalizacji kwadratowej \eqref{w_zad_opt_npltp}. Iteracje wewnętrzne zostają przerwane jeżeli różnica między przyszłymi przyrostami sygnałów sterujących obliczonymi w dwóch kolejnych iteracjach wewnętrznych jest mała, to znaczy gdy \begin{equation} \left\|\triangle\boldsymbol{u}^t(k)-\triangle\boldsymbol{u}^{t-1}(k)\right\|^2<\delta_{\mathrm{u}}\nonumber \end{equation} przy czym wielkość \(\delta_{\mathrm{u}}>0\) jest dobierana eksperymentalnie. Pierwszy element wektora zmiennych decyzyjnych \(\triangle\boldsymbol{u}^t(k)\) jest stosowany do sterowania. W kolejnej iteracji, po aktualizacji pomiarów, horyzont zostaje przesunięty jeden krok do przodu i cała procedura obliczeniowa zostaje powtórzona.
Aby ograniczyć nakład obliczeń można zastosować suboptymalne algorytmy regulacji predykcyjnej w wersji analitycznej, w których nie ma potrzeby optymalizacji kwadratowej [T07, T02]. Pomijając ograniczenia, zadanie optymalizacji algorytmu NPL \eqref{w_zad_opt_sl} można sprowadzić do problemu optymalizacji funkcji kwadratowej [Ł14, Ł09a] \begin{equation} \min\limits_{\triangle \boldsymbol{u}(k)} \left\{J(k)=\left\| \boldsymbol{y} ^{\mathrm{zad}}(k)-\boldsymbol{G}(k)\triangle\boldsymbol{u}(k)-\boldsymbol{y}^{0}(k)\right\|^{2}_{\boldsymbol{M}} + \left\| \triangle\boldsymbol{u}(k)\right\| ^{2}_{\boldsymbol{\Lambda}}\right\} \nonumber %\label{w_zad_optnpl_anal} \end{equation} Gradient powyższej funkcji kwadratowej względem zmiennych decyzyjnych algorytmu ma postać \begin{equation}\frac{\mathrm{d}J(k)}{\mathrm{d}\triangle \boldsymbol{u}(k)}=-2\boldsymbol{G}^{\mathrm{T}}(k)\boldsymbol{M}(\boldsymbol{y} ^{\mathrm{zad}}(k)-\boldsymbol{y}^{0}(k))+2\boldsymbol{\Lambda}\triangle \boldsymbol{u}(k)\end{equation} W punkcie optymalnym gradient zeruje się, z czego łatwo otrzymuje się wektor przyszłych przyrostów sygnałów sterujących \begin{equation} \triangle \boldsymbol{u}(k)=\boldsymbol{K}(k)(\boldsymbol{y} ^{\mathrm{zad}}(k)-\boldsymbol{y}^{0}(k))\label{w_du_K_yzady0_nplanal} \end{equation} gdzie \begin{equation} \boldsymbol{K}(k)=(\boldsymbol{G}^{\mathrm{T}}(k)\boldsymbol{M}\boldsymbol{G}(k)+\boldsymbol{\Lambda})^{-1}\boldsymbol{G}^{\mathrm{T}}(k)\boldsymbol{M} \nonumber %\label{w_K_nplanal} \end{equation} jest macierzą o wymiarach \(N_{\mathrm{u}}\times N\).
Walgorytmie analitycznym nie ma potrzeby obliczania całego wektora przyszłych przyrostów sygnałów sterujących \(\triangle \boldsymbol{u}(k)\), a jedynie \(n_{\mathrm{u}}\) jego pierwszych elementów, przyrostów sygnałów sterujących dla aktualnej iteracji. Zamiast wzoru (\ref{w_du_K_yzady0_nplanal}) stosuje się zależność \begin{equation} \triangle u(k|k)=\boldsymbol{K}^{n_{\mathrm{u}}}(k)(\boldsymbol{y} ^{\mathrm{zad}}(k)-\boldsymbol{y}^{0}(k))\label{w_du_Knu_yzady0_nplanal} \end{equation} gdzie macierz \(\boldsymbol{K}_{n_{\mathrm{u}}}(k)\) składa się z \(n_{\mathrm{u}}\) pierwszych wierszy macierzy \(\boldsymbol{K}(k)\) (dla procesu o jednym wejściu i jednym wyjściu \(n_{\mathrm{u}}=1\)). Macierz ta jest obliczana na bieżąco w każdej iteracji algorytmu, gdyż zależy od macierzy dynamicznej \(\boldsymbol{G}(k)\), zawierającej współczynniki odpowiedzi skokowej modelu zlinearyzowanego w aktualnym punkcie pracy. Do obliczania macierzy \(\boldsymbol{K}_{n_{\mathrm{u}}}(k)\) stosuje się wybraną numerycznie poprawną metodę wyznaczania macierzy odwrotnej, np. rozkład (faktoryzację) LU (ang. Low-Upper) macierzy, a następnie rozwiązuje się prosty układ równań liniowych [Ł09a].
Obliczone rozwiązanie może nie spełniać istniejących ograniczeń wartości oraz maksymalnej szybkości zmian sygnału sterującego \begin{equation} u^{\min}\le u(k|k) \le u^{\max}, \ -\triangle u^{\max}\le \triangle u(k|k) \le \triangle u^{\max}\nonumber \end{equation} Dlatego też obliczony wektor \(\triangle u(k|k)\) zostaje rzutowany na zbiór rozwiązań dopuszczalnych, wyznaczony przez ograniczenia. Procedura rzutowania jest następująca \begin{align} &\mathrm{jeżeli} \ \triangle u(k|k)<-\triangle u^{\max}, \ \mathrm{przyjąć} \ \triangle u(k|k)=-\triangle u^{\max}\nonumber\\ &\mathrm{jeżeli} \ \triangle u(k|k)>\triangle u^{\max}, \ \mathrm{przyjąć} \ \triangle u(k|k)=\triangle u^{\max}\nonumber\\ &\mathrm{obliczyć} \ u(k|k)=\triangle u(k|k)+u(k-1)\nonumber\\ &\mathrm{jeżeli} \ u(k|k)< u^{\min}, \ \mathrm{przyjąć} \ u(k|k)= u^{\min}\nonumber\\ &\mathrm{jeżeli} \ u(k|k)> u^{\max}, \ \mathrm{przyjąć} \ u(k|k)= u_n^{\max}\nonumber \end{align} Ograniczenia dotyczą tylko sygnałów sterujących dla aktualnej iteracji, gdyż w przeciwieństwie do algorytmów numerycznych nie oblicza się sekwencji przyrostów sygnałów sterujących na całym horyzoncie sterowania. Dla uproszczenia nie uwzględnia się również ograniczeń prognozowanych sygnałów wyjściowych.
Struktura analitycznego algorytmu NPL została przedstawiona na rys. 6, jest ona bardzo podobna do pokazanej na rys. 2 struktury klasycznego algorytmu NPL. Linearyzacja modelu neuronowego w aktualnym punkcie pracy i obliczenie nieliniowej trajektorii swobodnej odbywa się tak samo w obu algorytmach. W algorytmie analitycznym nie ma potrzeby optymalizacji kwadratowej, definiująca prawo regulacji \eqref{w_du_Knu_yzady0_nplanal} macierz \(\boldsymbol{K}_{n_{\mathrm{u}}}(k)\) zostaje obliczona w wyniku rozkładu macierzy.
Oczywiście, zamiast analitycznego algorytmu NPL możliwe jest również zastosowanie analitycznego algorytmu SL, jedyna różnica polega na użyciu do obliczania trajektorii swobodnej zlinearyzowanego, a nie nieliniowego, modelu procesu. Można również wyprowadzić algorytmy analityczne z bardziej złożoną linearyzacją trajektorii (algorytmy NPLT oraz NPLTP) [Ł14], zamiast prostej linearyzacji modelu w punkcie pracy. Pomijając ograniczenia, kwadratowe zadanie optymalizacji algorytmu NPLT \eqref{w_zad_opt_nplt} można zapisać w postaci \begin{align} &\begin{array}{l} \min\limits_{\triangle\boldsymbol{u}(k)} \Big \{ J(k)=\big\| \boldsymbol{y} ^{\mathrm{zad}}(k)-\boldsymbol{H}(k)\boldsymbol{J}\triangle \boldsymbol{u}(k)-\hat{\boldsymbol{y}}^{\mathrm{traj}}(k)\end{array}\nonumber\\ &\qquad \qquad \qquad \ \ -\boldsymbol{H}(k)(\boldsymbol{u}(k-1)-\boldsymbol{u}^{\mathrm{traj}}(k))\big\|^{2}_{\boldsymbol{M}}+ \left\|\triangle\boldsymbol{u}(k)\right\| ^{2}_{\boldsymbol{\Lambda}} \Big\} %\label{w_zad_optnplt_anal} \nonumber \end{align} Z warunku zerowania się gradientu minimalizowanej funkcji, otrzymuje się optymalny wektor przyrostów sygnałów sterujących \begin{equation} \triangle \boldsymbol{u}(k)=\boldsymbol{K}(k)(\boldsymbol{y} ^{\mathrm{zad}}(k)-\hat{\boldsymbol{y}}^{\mathrm{traj}}(k))+\boldsymbol{K}_{\mathrm{u}}(k)(\boldsymbol{u}(k-1)-\boldsymbol{u}^{\mathrm{traj}}(k)) \nonumber \end{equation} gdzie \begin{align} \boldsymbol{K}(k)&=(\boldsymbol{J}^{\mathrm{T}}\boldsymbol{H}^{\mathrm{T}}(k)\boldsymbol{M}\boldsymbol{H}(k)\boldsymbol{J}+\boldsymbol{\Lambda})^{-1}\boldsymbol{J}^{\mathrm{T}}\boldsymbol{H}^{\mathrm{T}}(k)\boldsymbol{M}\nonumber\\ \boldsymbol{K}_{\mathrm{u}}(k)&=-\boldsymbol{K}(k)\boldsymbol{H}(k) %\label{w_K_Ku_npltanal} \nonumber \end{align} są macierzami o wymiarach \(N_{\mathrm{u}}\times N\) oraz \(N_{\mathrm{u}}\times N_{\mathrm{u}}\). Nie ma oczywiście potrzeby obliczania całego wektora przyszłych przyrostów sygnałów sterujących \(\triangle \boldsymbol{u}(k)\), a jedynie \(n_{\mathrm{u}}\) jego pierwszych elementów - przyrostów dla aktualnej iteracji \begin{equation} \triangle u(k|k)=\boldsymbol{K}^{n_{\mathrm{u}}}(k)(\boldsymbol{y} ^{\mathrm{zad}}(k)-\hat{\boldsymbol{y}}^{\mathrm{traj}}(k))+\boldsymbol{K}_{\mathrm{u}}^{n_{\mathrm{u}}}(k)(\boldsymbol{u}(k-1)-\boldsymbol{u}^{\mathrm{traj}}(k)) \nonumber \end{equation} gdzie macierz \(\boldsymbol{K}^{n_{\mathrm{u}}}(k)\) zawiera \(n_{\mathrm{u}}\) pierwszych wierszy macierzy \(\boldsymbol{K}(k)\), natomiast macierz \(\boldsymbol{K}^{n_{\mathrm{u}}}_{\mathrm{u}}(k)\) zawiera \(n_{\mathrm{u}}\) pierwszych wierszy macierzy \(\boldsymbol{K}_{\mathrm{u}}(k)\).
Poprzednia strona: 1. Zasada regulacji predykcyjnej | Następna strona: 3. Sieci neuronowe w regulacji predykcyjnej |