Materiały pomocnicze do ćwiczeń
z  Nowoczesnych metod sterowania
część 4: Sterowanie adaptacyjne i predykcyjne.
ĆW.9. Identyfikacja obiektów sterowania w trybie nadążnym.
ĆW.10. Adaptacyjne metody sterowania.
ĆW.11. Sterowanie predykcyjne.
Jan Leszczyński
Katedra Automatyki i Elektroniki
Politechnika Białostocka.
BIAA
YSTOK, 2001
2
ĆW.9. Identyfikacja obiektów sterowania w trybie nadążnym.
Porównywane są wybrane wersje sieci działań identyfikacji w trybie nadążnym.
W trakcie doświadczeń symulacyjnych badane są właściwości algorytmów rekurencyj-
nych identyfikacji na tle bieżącej estymacji parametrów zmieniającego się modelu ob-
serwowanego procesu dynamicznego, m.in. badana jest zdolność śledzenia tych
algorytmów za zmianami parametrów w modelu o określonej strukturze.
Istotną składową ćwiczenia jest także zasygnalizowanie problemów komunikatywnego
przedstawiania w czasie rzeczywistym dyspozytorowi procesu informacji uzyskiwanej
 na bieżąco przez cyfrowy system pomiarowo - sterujący.
PRZED ĆWICZENIEM:
Y
Y Należy pogłębić i uzupełnić wiedzę na temat łącznego modelowania, w konwencji dys-
kretnej, obiektu dynamicznego i oddziałujących nań zakłóceń. Powtarzając materiał z za-
kresu rekurencyjnych algorytmów identyfikacji, kładziemy nacisk na problem
dopasowania stosowanego algorytmu do postaci modelu procesu dynamicznego. Podsta-
wowe informacje na ten temat ujęto w zestawieniu A.
Pytania kontrolne:
Jaką postać przyjmuje rekurencyjna wersja estymatora MNK w zastosowaniu do identyfi-
kacji modelu dynamicznego w strukturze ARX?
W jaki sposób wprowadza się w tym algorytmie preferencje dla informacji najbardziej
aktualnej?
Jakie dodatkowe zadania obliczeniowe występują przy opisie dynamiki obiektu i zakłóceń
przy zastosowaniu modelu ARMAX? Przedstaw program obliczeń jednego etapu realiza-
cji tego algorytmu.
Czy istnieje rekurencyjna wersja estymatora zmiennej instrumentalnej? W jaki sposób
obliczenia iteracyjne w tej wersji są ujmowane w ramy uwarunkowań czasowych ?
3
Y
Y Należy przypomnieć podstawowe informacje na temat pakietu Identyfikacji Ma-
tlaba. Przedstawiono je zwięz
le w zestawieniu D w dalszej części instrukcji. Wskaza-
nym jest także zapoznanie się z tekstem pomocy do rekurencyjnych m- funkcji
specjalizowanego pakietu Identyfikacji( np.  ) w środowisku
programowym Matlaba- w nazwach tych funkcji, rozpoczynających się literą  r wid-
nieją rodzaje stosowanych struktur przy modelowaniu procesów dynamicznych. W dalszej
części instrukcji znajdują się m.in. programowane sekwencje działań umożliwiających
przeprowadzenie symulacji eksperymentu ukierunkowanego na identyfikację modelu
obiektu dynamicznego przy zastosowaniu obliczeń rekurencyjnych.
Pytania kontrolne:
W jaki sposób można wprowadzić m- funkcjach i warunki po-
czątkowe dotyczące macierzy  P i wektora estymowanych parametrów?
W jakich wersjach, w rekurencyjnych algorytmach, można zadeklarować celowość zasto-
sowania czynnika związanego z zapominaniem zdezaktualizowanych danych?
METODYKA BADAC
SYMULACYJNYCH I PROGRAMOWANIE:
Identyfikowany obiekt jest zadawany za pośrednictwem modelu analogowego i repre-
zentowany w obliczeniach symulacyjnych przy pomocy modelu dyskretnego w postaci linio-
wego równania różnicowego wiążącego wielkość wyjściową z oddziaływaniem sterującym
i zakłócającym. Przy zastosowaniu modelu w strukturze ARX wyrażenie na wyznaczenie
bieżącego poziomu wyjścia przyjmuje w związku z tym postać:
w której:
A i B są wielomianami (stopnia odpowiednio na, nb) licznika i mianownika transmitan-
cji dyskretnej obiektu,
yw, uw reprezentują aktualne, czasowo uporządkowane, odpowiednio na i (nb+1) wy-
miarowe, wektory wyjść i sterowań.
Przez pzoznaczono poziom zakłóceń.
Czynnik zakłócający jest tu wyrażany przy pomocy aktualnej realizacji gaussowskiej zmien-
nej losowej mnożonej przez pz. Przy zastosowaniu modelu ARMAX sytuacja nieco się
4
komplikuje; czynnik ten reprezentowany jest iloczynem skalarnym dwu wektorów. Jeden z
nich to wektor współczynników wielomianu C, drugim natomiast sekwencja ostatnich, staty-
stycznie niezależnych losowań oddziaływania wejściowego w torze zakłóceń. Jest ono natu-
ry czysto probabilistycznej  przebieg reprezentujący biały szum jest generowany za
pośrednictwem standardowego generatora liczb losowych   .
Oddziaływanie testujące od strony wejścia sterowanego obiektu ma natomiast z założe-
nia charakter pseudo-stochastycznego sygnału PRBS (ang. Pseudo Random Binary Signal ) ,
np. [1]. Program generatora PRBS w środowisku programowym Matlaba może przyjąć
następującą postać:
sbp = rej(:,1);
Badanie właściwości algorytmów identyfikacji obiektów dynamicznych możemy
przeprowadzić najsprawniej posługując się gotowymi składnikami, pakietu oprogramowania
Identyfikacji (zestawienie D). Prostota programu doświadczenia symulacyjnego wyni-
ka z możliwości zastosowania specjalizowanej m-funkcji  idsim . Funkcja ta wymaga
uprzednio wygenerowania modelu w standardowej ( w ramach danego pakietu) postaci
 tetha . Przytoczymy charakterystyczny przykład programu symulacji procesu badania
obiektu dynamicznego opisywanego modelem ARMAX  zakłada się, że widoczne tu para-
metry  pst ,  pz zostały ustalone przed jego wywołaniem:
u= pst * sign(randn(100,1)); % symetryczny sygna" testuj" cy
e= pz*randn(100,1); % generator zak" óce" amplitudy  pz
th0= poly2th([1 -1.5 0.7],[0 1 0.5],[1 -1 0.2]); % model ARMAX
5
y = idsim([u e],th0); % instrukcja proces symulacji
z = [y u]; idplot(z), % wizualizacja sygna" ów
Przedstawimy również typową sekwencję rekurencyjnych obliczeń estymacji parametrów
modelu w strukturze najmniejszych kwadratów (ARX).
% wst" pny etap oblicze" z przygotowaniem do wizualizacji wyników
[th,yh,p,phi] = rarx(z(1,:),[2 2 1],'ff',0.98);
plot(1,th(1),'*',1,th(2),'+',1,th(3),'o',1,th(4),'*'),
axis([1 100 -2 2]), hold on; % wymiary na osiach wspó" rz" dnych
% obliczenia rekurencyjne wraz z wizualizacj" wyników
for k = 2:100
[th,yh,p,phi] = rarx(z(k,:),[2 2 1],'ff',0.98,th',p,phi);
plot(k,th(1),'*',k,th(2),'+',k,th(3),'o',k,th(4),'*')
end
Znaczenie argumentów wejściowych i wyjściowych stosowanych m-funkcji przedstawio-
ny jest w ich  helpach . Celowym jest także zapoznanie się z całością programu ilustrującego
działanie algorytmów rekurencyjnych   , z którego pochodzi dana sekwencja in-
strukcji. Podkreślmy, że program identyfikacji przeprowadzany jest tu w trybie off line. Za-
kłada się przed jego rozpoczęciem obecność w pamięci komputera otrzymanego w trakcie
pomiarów lub symulacji ciągu czasowego wyjść i sterowań obiektu: z = [y,u].
6
PROGRAM ĆW. 9
1. Przeprowadzamy doświadczenie symulujące proces identyfikacji obiektu dynamicznego
jednym z algorytmów rekurencyjnych przy następujących założeniach:
W trakcie doświadczenia symulacyjnego stosowany jest sygnał testujący w postaci pseu-
do-przypadkowego sygnału PRBS ustalanej długości i amplitudy. y
ródłem zakłóceń jest
natomiast sygnał białego szumu o rozkładzie normalnym.
Dane o symulowanych przebiegach opracowywane są przy pomocy znajdującej się w pa-
kiecie tematycznym Identyfikacji, m- funkcji  lub  Należy
zapamiętać przebiegi zarówno wielkości wejściowej u(k) i wyjściowej y(k) badanego
obiektu (w postaci macierzy z=[y u]) jak i przebiegów estymowanych w trybie nadąż-
nym parametrów jego modelu cest(k).
Rezultatem doświadczenia symulującego proces identyfikacji obiektu są (m.in. wyrażone
graficzne) przebiegi estymacji parametrów jego modelu. Początkowe wartości oszacowań
parametrów we wszystkich doświadczeniach przyjmujemy jako zerowe.
Model badanego procesu dynamicznego jest reprezentowany transmitancją Laplace a
pierwszego lub drugiego rzędu z opóz
nieniem. Przykłady zalecanych struktur modelu
i danych liczbowych:
�
�
� � �� �
�
� � �� �
Np.: � � �
7
 Tor zakłóceń w modelu reprezentuje z założenia niezależnie definiowany wielomian
C(z-1), np. C(z-1)=1; (otrzymujemy wówczas model ARX), C(z-1)=1+0,5z-1;
C(z-1)= 1+0,7z-1+0,1z-1; ... itp.
2. Należy zaplanować i przeprowadzić doświadczenie symulacyjne z procesem identyfikacji
parametrów modelu w celu ustalenia zależności między szybkością zbieżności algorytmów
rekurencyjnych identyfikacji, a stosunkiem mocy sygnału testującego do wariancji zakłó-
ceń zniekształcających pomiary wielkości wyjściowej badanego obiektu.
Wskazówka. Szybkość zbieżności algorytmu możemy mierzyć np. przy pomocy liczby taktów
procesu przejściowego wymaganych do momentu umownego ustalenia się oszacowań pa-
rametrów modelu. Program badań powinien umożliwić między innymi udzielenie miaro-
dajnej odpowiedzi na pytanie: Jaki stosunek sygnału testującego do zakłóceń można
w sytuacji zdeterminowanej założeniami danego eksperymentu symulacyjnego uznać za
minimalnie - dopuszczalny ?
3. Należy zmodyfikować doświadczenie symulacyjne w ten sposób by można było zaobser-
wować jakościowo proces śledzenia estymatora rekurencyjnego za wolno zmieniającym
się parametrem niestacjonarnego modelu.
Wskazówka. W tym celu zakłada się 30% zmianę o charakterze liniowym jednego lub dwu
parametrów transmitancji ciągłej badanego procesu ( np. współczynnika tłumienia, lub
współczynnika wzmocnienia) w ciągu 200 taktów pomiarowych. W stosowanych algo-
rytmach rekurencyjnych wprowadzamy mechanizm zapominania, deklarując trzeci argu-
ment wejściowy w postaci: (ang. forgetting factor ) i wprowadzając (w czwartym
argumencie) wybraną wartość współczynnika zapominania. Należy zbadać, definiując sa-
modzielnie określony wskaz
nik jakości, jaki wpływ na rezultaty obliczeń wywierają zmia-
ny współczynnika zapominania w granicach od  do  . Jaki
współczynnik zapominania należałoby uznać w danej sytuacji za optymalny?
4. Na podstawie znanego równania obiektu i jednego z równań regulatorów otrzymanych
w poprzednich badaniach zamodelować układ stabilizujący wartość regulowaną na okre-
ślonym poziomie, a następnie powtórzyć proces identyfikacji prowadząc obserwacje
w układzie zamkniętym, przy tych samych warunkach początkowych, zakłóceniach
8
i addytywnym sygnale testującym co w punkcie 1 ćwiczenia. Czy różnice przebiegów es-
tymacji parametrów modelu obiektu są istotne? Jakie są główne przyczyny ewentualnych
rozbieżności?
ZADANIE NADOBOWI
ZKOWE jest związane z przekazywaniem informacji w trybie
nadążnym. Zaleca się zmodyfikować program doświadczenia symulacyjnego w ten sposób
by można było zaobserwować jakościowo proces estymacji wolno zmieniających się pa-
rametrów niestacjonarnego modelu obiektu sterowania. W ogólniejszym sformułowaniu
jest to zadanie bliskie zaprogramowaniu i przetestowanie programu, którego celem jest
zilustrowanie procesu przekazywania metodą graficzną dyspozytorowi procesu  na bieżą-
co informacji o aktualnych wynikach identyfikacji i jakości aktualnego modelu obserwo-
wanego procesu.
Wskazówka. Jednym ze sposobów rozwiązania tego zadania jest uzyskanie przesuwającego
się w czasie przebiegu ocen parametrów modelu w określonym ( np. kilkanaście ostatnich
ocen) przedziale czasowym. Przedział ten powinien przesuwać się wraz z upływem czasu
w ten sposób, by chwila bieżąca zajmowała ustalone, zawsze to samo miejsce. Nazwa an-
gielska wykresu tego typu: scrolling chart.
Celowym jest również umieszczenie na innym nadążnie aktualizowanym wykresie
przebiegu wielkości mierzonej y(k) i jej oceny yest(k) uzyskiwanej na podstawie
aktualnego modelu procesu. W jakiej sytuacji celowym mogłoby być także przedstawianie
w trybie nadążnym dyspozytorowi procesu równolegle kilku wykresów z różnymi skalami
czasowymi?
SPRAWOZDANIE POWINNO ZAWIERAĆ:
Programy przeprowadzonych doświadczeń symulacyjnych i graficzną prezentację ilu-
stratywnych przebiegów ocen estymowanych parametrów, np. wskazanym jest porówna-
nie na jednym wykresie nadążnych wyników estymacji dla dwu skrajnych, stosowanych
w trakcie badań współczynników zapominania.
9
 Komentarze i wnioski z przeprowadzonych w ćwiczeniu badań. Czy stawiane cele zo-
stały osiągnięte? Przedstaw uzasadnioną wynikami doświadczeń symulacyjnych odpo-
wiedz
na pytania stawiane w punkcie drugim i trzecim programu ćwiczenia.
10
ĆW 10. Adaptacyjne metody sterowania.
W trakcie obliczeń i doświadczeń symulacyjnych prowadzonych w środowisku pro-
gramowym Matlaba/Simulink badane są podstawowe wersje mechanizmów re-
gulacji adaptacyjnej: układy adaptacji pośredniej z identyfikacją obiektu sterowania
i układy z wbudowanym mechanizmem bezpośredniej identyfikacji parametrów regula-
tora. Uwypuklana jest przy tym problematyka nadążania zmian własności regulatora
za zmianami parametrów sterowanego obiektu i ocena czasów realizacji poszczegól-
nych segmentów oprogramowania w okresowo powtarzanym cyklu obliczeń.
PRZED ĆWICZENIEM:
Y
Y Należy powtórzyć korzystając z notatek z wykładu i literatury (np. [2]) i kierując się za-
wartością zestawienia B podstawowe elementy teorii układów adaptacyjnych.
Pytania kontrolne:
Jakie różnice występują w schematach blokowych układów regulacji adaptacyjnej
w wersji pośredniej i bezpośredniej. Porównaj na tym tle cechy charakterystyczne tych
układów.
Wymień główne składowe sieci działań realizujących wymienione mechanizmy sterowania
adaptacyjnego w przypadku, gdy celem jest ustalenie biegunów układu regulacji.
W jaki sposób można ocenić jakość układu regulacji adaptacyjnej?
Badanie właściwości algorytmów rekurencyjnych regulacji adaptacyjnej
przeprowadzamy programując cyklicznie powtarzaną czteroetapową sekwencję obliczeń
obejmującą:
A. symulację pomiaru wielkości wyjściowej obiektu sterowania; jest ona generowana na
podstawie znanych aktualnych wektorów oddziaływania sterującego i zakłócającego;
B. aktualizację wyniku estymacji parametrów modelu sterowanego obiektu jednym z
algorytmów identyfikacji rekurencyjnej (przy stosowaniu mechanizmu adaptacji bezpo-
średniej ta sekwencja obliczeń jest opuszczana);
11
C. wyznaczenie i wprowadzenie nowych, skorygowanych wartości parametrów stosowa-
nego regulatora;
D. wyznaczenie nowego poziomu sterowania w wyniku zrealizowania obliczeń wynikają-
cych z zaktualizowanego równania regulatora  oddziaływanie to jest miarodajne w na-
stępnym takcie pracy układu.
Wyniki działania każdego z wymienionych etapów obliczeń są zachowywane w celu ich
przekazania po zakończeniu doświadczenia symulacyjnego. Przedstawimy kilka komentarzy
do wymienionych etapów obliczeń.
A. Strukturę modelu obiektu sterowania wraz z torem zakłóceń stosujemy analogiczną
do badanej w poprzednim ćwiczeniu; informację wyjściową o modelu w postaci analogowej
należy przekształcić, przy uwzględnieniu zakłóceń, do równania różnicowego typu ARX lub
ARMAX. Przed rozpoczęciem procesu symulacji należy ustalić wartości startowe wektorów
yw, uw ( i ewentualnie epsw) widocznych w tym równaniu:
Zwykle przyjmuje się je jako zerowe. Zakłócenia przypadkowe wprowadzane są na bieżąco
przy pomocy generatora liczb losowych Matlaba, a oddziaływania sterujące może być uzu-
pełniane addytywnym sygnałem testującym, np. typu PRBS z generatora binarnego.
W trakcie ćwiczenia może być wykorzystywana zarówno stacjonarna jak i niestacjonarna
wersja obiektu sterowania.
B. Algorytm identyfikacji wybieramy na podstawie wyników badań w ćwiczeniu po-
święconym algorytmom rekurencyjnym identyfikacji. Istotnym elementem jest tu m.in. racjo-
nalny wybór  współczynnika zapominania . Istnieje przy tym konieczność wyboru
i wprowadzenia do przestrzeni roboczej Matlaba określającej kowariancję ocen, startowej
macierzy P(0)i startowego wektora estymowanych parametrów w modelu np. przy pomocy
instrukcji:
C. Dostrajanie parametrów regulatora może odbywać się metodą najmniejszych
kwadratów, w metodyce bezpośredniej, lub poprzez ich cykliczne wyliczanie po korekcie
wyników identyfikacji sterowanego obiektu. Z równania różnicowego regulatora:
12
wynika, że koniecznym jest wyliczanie w trybie nadążnym (nr+ns+nt+2) współczynników
wielomianów R, S, T. Symulując mechanizm adaptacji pośredniej stosujemy jedną z me-
tod syntezy regulatora ujętą w postaci programu m-funkcji opracowanej w ramach poprzed-
nich ćwiczeń na temat regulatorów cyfrowych (PPC, PZPC, lub MVC). W układach regulacji
adaptacyjnej z bezpośrednią identyfikacją parametrów regulatora stosowany jest natomiast
mechanizm rekurencji analogiczny do omawianego przy algorytmach identyfikacji obiektu.
Proces adaptacji rozpoczyna się od nie nastrojonego regulatora w którym zakładamy po-
czątkowe wektory parametrów R, S, T rozmiaru wynikającego z przyjętych a priori stopni
wielomianów, np. T=1; S=[.01 .01]; R = [1 -.01 .01]. Przy ustalaniu warun-
ków początkowych należy pamiętać o tym, by wektory początkowe z równania
regulatora nie były sprzeczne z wektorami początkowymi przyjętymi dla obiektu sterowania.
D. Sygnał sterujący jest wyznaczany z równania regulatora. Przypomnijmy, że wyrażany
jest on za pośrednictwem trzech wielomianów regulatora: S(z-1), T(z-1), R(z-1) okre-
ślających związek między wielkościami: zadaną w(k), regulowaną y(k) i sterującą
u(k):
R(z-1) u(k) = T(z-1) w(k) - S(z-1) y(k).
Do równania tego wprowadza się najbardziej aktualne, wyznaczone wcześniej w tej sa-
mej sekwencji działań, wartości jego parametrów. Sygnał wartości zadanej w(k)generujemy
w postaci impulsów prostokątnych o okresie kilkudziesięciu taktów ponieważ w przebiegach
przejściowych zawarta jest podstawowa informacja o zmieniających się z założenia własno-
ściach dynamicznych układu regulacji.
Podstawowy sygnał sterujący może być potencjalnie uzupełniany sygnałem testującym,
np. z generatora binarnego PRBS, lub poprzez wprowadzenie do programu instrukcji:
Każde wyliczenie nowego sterowania implikuje zmianę wektorów uw i uv w równa-
niach obiektu i regulatora, a w następstwie zmianę wektorów yw i yv. W programie sy-
mulacyjnym ujmuje się to w odpowiednich instrukcjach przesunięć.
13
Wskaz
niki oceny jakości badanego układu regulacji i wizualizacja przebiegów proce-
sów sterowania uzyskiwane są po zakończeniu obliczeń symulacyjnych. Mogą być one zo-
rientowane zarówno tradycyjnie, np. w powiązaniu z oceną procesów przejściowych,
wynikających przy skokowych przestawieniach wartości zadanej lub zdolności wytłumiania
zakłóceń, jak i mniej standartowo.
Jedna z wersji oceny niestandardowej może np. polegać na  rozłącznej ocenie prze-
biegu procesów identyfikacji i dostrajania regulatora przy wykorzystaniu znajomości a priori
rzeczywistych parametrów obiektu i dopasowanych, wyliczonych off line zgodnie ze stoso-
waną metodą syntezy, parametrów regulatora. Przykłady oprogramowania obliczeń tego typu
pokazano (m.in. komentarze w tekście programu) w programach badania układów regulacji
adaptacyjnej w zestawieniach D w danej instrukcji.
PROGRAM ĆW. 10
1. Powtórzyć doświadczenie symulujące proces identyfikacji obiektu dynamicznego jednym
z algorytmów rekurencyjnych z poprzedniego ćwiczenia. W trakcie doświadczenia odno-
towujemy w pamięci komputera przebiegi czasowe estymowanych parametrów obiektu.
2. Sięgając po program syntezy regulatora PPC lokalizującego bieguny wynikowego układu
regulacji stosowany w jednym z poprzednich ćwiczeń wyznaczamy wielomiany z równa-
nia regulatora: T(z-1)w(k)  S(z-1)y(k) = R(z-1)u(k), odpowiadające charak-
terystycznym punktom przebiegów czasowych estymowanych parametrów z poprzedniego
punktu, a następnie wyznaczamy prawdziwe położenie biegunów układów regulacji przy
uzyskanych równaniach regulatora. Odnotowujemy graficznie trajektorie przesuwania się
biegunów. Po ilu taktach obliczeń własności układu regulacji oceniane przez pryzmat po-
łożenia biegunów stabilizują się?
3. Należy opracować i przetestować główny składnik programu symulujący proces dostraja-
nia parametrów regulatora do nieznanego a priori, identyfikowanego modelu sterowanego
obiektu. W programie tym należy umieścić w ramach jednej z instrukcji cyklu ( for ,
 while , itp.) sekwencję obliczeń obejmującą:
Etap1. Obliczenie aktualnej wartości wielkości regulowanej na podstawie modelu obiektu:
14
Ay(k)=z^(-d)Bu(k)+C*e(k).
Etap2. Uaktualnienie wyniku identyfikacji przy zastosowaniu odpowiedniej m- funkcji Ma-
tlaba(np.  rarx ,  itp.).
Etap3. Wyznaczenie nowych parametrów regulatora przy pomocy jednej z przeznaczonych
do tego celu, znanej z poprzednich ćwiczeń m- funkcji (np. ppc ,  pzp1 ,  mvc ).
Etap4. Obliczenie z aktualnego równania regulatora: Tw(k)  Sy(k) = Ru(k); war-
tości sterowania obowiązującej w następnym cyklu obliczeń adaptacji rozpoczyna się od
nie nastrojonego regulatora i uzupełnienie jej aktualną wartością pomocniczego oddziały-
wania z generatora sygnałów typu PRBS.
W trakcie testów programu odnotowujemy czasy realizacji poszczególnych składowych
obliczeń i rejestrujemy kilka przebiegów dostrajania, reprezentatywnych dla różnych po-
ziomów zakłóceń. W rezultacie badań symulacyjnych należy sformułować wypowiedz
na
następujące tematy:
Czy czas trwania procesu przejściowego dostrajania regulatora zależy od poziomu za-
kłóceń?
Czy obecność dodatkowego sygnału testującego jest konieczna?
Wskazówka. Możemy skorzystać przy tym z przykładu sekwencji działań tego typu umiesz-
czonej w ramach programu  nms_10 , przytoczonego w zestawieniu D. W programie tym
skompletowano program doświadczenia symulującego proces dostrajania parametrów re-
gulatora wprowadzając m.in. generator wartości zadanej w postaci ciągu impulsów prosto-
kątnych o okresie kilkudziesięciu taktów i końcowe procedury obliczenia wskaz
ników
jakościowych przebiegów.
SPRAWOZDANIE POWINNO ZAWIERAĆ teksty opracowanych w trakcie realizacji
ćwiczenia programów, reprezentatywne przebiegi kluczowych wielkości w trakcie obli-
czeń i inne istotniejsze wyniki przeprowadzonych badań symulacyjnych, a także zwięzłe
wypowiedzi na tematy poruszane w przedstawionych kursywą pytaniach.
15
ĆW 10. Sterowanie predykcyjne.
Badany jest, w trakcie doświadczeń symulacyjnych, wpływ przyjmowanych w procesie
projektowania założeń wyjściowych i podstawowych parametrów na rezultaty syntezy
układu sterowania predykcyjnego. Program ćwiczenia przewiduje aktywne zapoznanie
się z zawartością specjalizowanego tematycznie pakietu programowego o nazwie: Mo-
del Predictive Control Toolbox poprzez prześledzenie typowych sekwencji
instrukcji spotykanych w programach reprezentatywnych dla jego zastosowań.
Porównane są także rezultaty kompensacji deterministycznego oddziaływania zakłó-
cającego w klasycznym układzie regulacji i w ramach sterowania predykcyjnego.
PRZED ĆWICZENIEM:
Y
Y Należy powtórzyć, kierując się zawartością zestawienia A w dalszej części opracowania,
materiał teoretyczny na temat syntezy sterowania predykcyjnego.
Pytania kontrolne:
W jaki sposób można przedstawić schemat układu sterowania predykcyjnego wyrażając
(blokowo) główne zadania obliczeniowe cyklicznie rozwiązywane przez system cyfrowy
w trakcie wypracowywania oddziaływania sterującego?
W jaki sposób wyznaczana jest zwykle pożądana trajektoria odniesienia?
Jaką rolę odgrywają pojęcia horyzontu predykcji i horyzontu sterowania.
Jaki kwadratowy wskaz
nik jakości stosowany jest zwykle do oceny jakości sterowania
predykcyjnego ?
Y
Y Należy zapoznać się z zawartością pakietu (tzw. toolboksu) Matlaba na temat stero-
wania predykcyjnego Model Predictive Control(w skrócie MPC). Jeżeli w
posiadana wersja Matlaba dysponuje tym pakietem oprogramowania wystarczy wywołać
w tym celu:  help mpccmds , a następnie  help mpccon ,  help mpcsim . Po-
winniśmy także przeanalizować gotowe przykłady projektowania zawartych np. w m-
funkcjach:  mpcconex ,  smpccoex ,  cmpcex . Informacje wprowadzające do te-
matu zawarto w zestawieniu E.
16
Pytania kontrolne:
Na jakie bloki funkcjonalne zostały podzielone m- funkcje danego pakietu?
Jakie typy modeli obiektów sterowania stosowane są tu standardowo?
Przedstaw argumenty wejściowe m- funkcji syntezy  mpccon regulatora predykcyjne-
go. W jaki sposób generowana jest tu trajektoria odniesienia?
Czy pakiet zawiera programy procesu symulacji układu sterowania predykcyjnego
i przedstawienia (m.in. wizualizacja) głównych jego wyników?
METODYKA BADAC
I PROGRAMOWANIE:
W trakcie badań symulacyjnych przeprowadzanych w ćwiczeniu posługujemy
się, podobnie jak poprzednio, przygotowanymi segmentami oprogramowania związanymi:
1. z wyrażeniem modelu obiektu sterowania w formacie dopasowanym do wymaganego
przez m- funkcję syntezy regulatora;
2. z programem syntezy regulatora lub realizacji innej metody wyznaczenia oddziaływania
sterującego;
3. z symulacją działania uzyskanego układu sterowania i obliczaniem wskaz
ników jako-
ściowych, a także wizualizacją typowych przebiegów przejściowych.
Dalej przedstawimy przykłady konstrukcji segmentów oprogramowania wymienionych ty-
pów. Opracowano je korzystając z gotowych m- funkcji pakietu: Model Predictive
Control(w skrócie mpccmds).
Konkretne zadanie badawcze realizujemy zmieniając nadrzędnie, przy pomocy jednej
z instrukcji cyklu wartość badanego parametru w interesującym nas przedziale (przykłady
zawarto w zestawieniu C). Na program ćwiczenia składa się realizacja kilku zadań badaw-
czych.
Model sterowanego obiektu należy przetransformować do formatu stosowanego
w ramach pakietu mpccmds. Zakładając, że jest to oparty na transmitancji model obiektu
SISO, analogiczny do stosowanego w poprzednich ćwiczeniach, przekształcenie to przepro-
wadzamy stosując np. dwie m- funkcje: najpierw:   , a następnie:   .
17
W międzyczasie wymagane jest (argumenty danych funkcji) ustalenie czasu próbkowania,
liczby wyjść w systemie, a także efektywnej długości odpowiedzi skokowej obiektu.
Przedstawimy przykład segmentu oprogramowania związanego z wypracowaniem modelu
sterowanego obiektu w odpowiednim formacie możemy (dane liczbowe pochodzą z wcze-
śniej stosowanego przykładu).
W rezultacie realizacji tego segmentu programu w przestrzeni roboczej Matlaba znaj-
dzie się model sterowanego obiektu w postaci dyskretnej charakterystyki skokowej wymaga-
nej w pakiecie MPC. Podkreślono w nim celowość zwrócenia uwagi na komunikat
automatycznie generowany przez m-funkcję:   . Komunikat ten informuje o osią-
ganej dokładności wyrażenia charakterystyki przy pomocy ograniczonego (przez zmienną
 tfinal ) ciągu jej wyrazów. W przypadku zbyt dużego procentu błędu koniecznym jest
przeprowadzenie korekty założeń i wydłużenie  tfinal .
Program syntezy regulatora stanowi następną fazę programowania. Przechodzimy
do niej po przyjęciu założenia o braku ograniczeń na wartość oddziaływania sterującego
i wielkości regulowanej w projektowanym układzie. Podstawę syntezy stanowi wówczas
standardowa, korzystająca z modelu uzyskiwanego w poprzedniej fazie obliczeń, m-funkcja
pakietu MPC o nazwie  mpccon . Przed wywołaniem  mpccon koniecznym jest ustalenie
parametrów sterowania predykcyjnego, tj. mierzonego liczbą taktów horyzontu predykcji P i
horyzontu sterowania M, a także wag ywt, uwt( uwaga: potencjalnie mogą to być wekto-
ry akcentujące problematykę następstw czasowych) ze wskaz
nika oceny jakości sterowania.
18
Uzyskiwanym rezultatem jest natomiast gotowy wektor wzmocnień składowych aktualnego
wektora sterowań. Segment oprogramowania związany z syntezą regulatora jest w związku z
tym bardzo prosty; obejmuje jedynie przygotowanie danych do wywołania i samo wywołanie
 mpccon :
Podkreślmy, że prostota rozwiązania zadania syntezy regulatora predykcyjnego w ramach
danego środowiska programowego może być złudną, gdyż sam wybór dostatecznie dobrych
w danym zadaniu parametrów projektowania nie zawsze jest sprawą prostą. Z drugiej strony
wybór ten jest istotnie ułatwiony w wyniku możliwości wielokrotnej symulacji skutków roz-
wiązywania zadania syntezy, przy zastosowaniu różnych wersji argumentów (
m- funkcji  mpccon .
Wyznaczanie przebiegów i wskaz
ników jakości sterowania stanowi od-
dzielny problem, również znacznie uproszczony w wyniku obecności gotowych, w ramach
omawianego pakietu oprogramowania, specjalizowanych m- funkcji analizy własności (m.in.
typowe przebiegi czasowe, charakterystyki częstotliwościowe i położenie biegunów) uzyski-
wanych układów sterowania. Szczególnie istotną sprawą jest tu umożliwienie bezpośredniej
symulacji działania układu regulacji przy zastosowaniu predykcyjnych metod syntezy od-
działywania sterującego, nawet jeżeli jest ono wyznaczane niezależnie w każdym takcie pracy
układu. Podkreślmy, że w tym trybie działania istnieje możliwość uwzględnienia niestacjo-
narności modelu sterowanego procesu, obecności ograniczeń sterowania, itp.
Przedstawimy przykład segmentu programowego oceny układu regulacji predykcyjnej.
W wyniku realizacji określonych tam instrukcji otrzymujemy przebiegi wielkości regulowa-
nej i oddziaływania sterującego w odpowiedzi na zadaną (tutaj skok jednostkowy) postać sy-
gnału wartości zadanej i ocenę jakości procesu przestawienia wartości zadanej wyrażoną
poprzez sumę kwadratów odchyłek w deklarowanym przedziale czasowym symulacji.
19
Dodajmy, że dysponując przebiegiem wielkości wyjściowej i oddziaływań sterujących
możemy łatwo wyznaczyć także inne wersje wskaz
nika jakości.
PROGRAM ĆW. 11:
1. Należy zaprogramować i zrealizować proces syntezy regulatora predykcyjnego MPC
(ang.  Model Predictive Controller) dla zadanej transmitancji obiektu sterowanego po-
sługując się odpowiednim pakietem oprogramowania narzędziowego Matlaba, a na-
stępnie przeprowadzić symulację uzyskanego układu regulacji.
Wskazówka. Pomocnym może być przy tym złożenie (połączone z odpowiednia modyfikacją)
przedstawionych wcześniej w instrukcji segmentów oprogramowania. Parametry syntezy
ustalamy wyjściowo zgodnie z wskazówkami literaturowymi i intuicyjnymi poglądami na
ten temat. Charakterystykę skokową sterowanego obiektu odtwarzamy z dokładnością po-
niżej 1%.
2. Należy ułożyć i zrealizować program badań symulacyjnych na temat wpływu długości
uwzględnianej w obliczeniach części charakterystyki skokowej (funkcji wagowej) na ja-
kość syntezy regulatora predykcyjnego. Jaka długość stosowanego ciągu czasowego jest
w badanej sytuacji wystarczająca? Jak ocenia się przy tym (w procentach) względną do-
kładność odtwarzanej charakterystyki?
3. Należy zaprogramować i zrealizować program badań symulacyjnych na temat doboru ho-
ryzontów predykcji P i sterowania M przy racjonalnym rozwiązaniu przydzielonego
(obiekt sterowany i postać wskaz
nika jakości układu regulacji) zadania regulacji predyk-
cyjnej. W wyniku badania powinien zostać przeprowadzony (uzasadnienie) wybór tych pa-
rametrów w rozwiązywanym zadaniu.
20
Wskazówka. W tym celu odnotowujemy wartość wskaz
nika uzyskanego przy wielokrotnym
rozwiązaniu, np. na podstawie programu z poprzedniego punktu umieszczonego w ramach
podwójnej instrukcji  for obejmującej elementy określonego zbioru {P, M}. Przy anali-
zie uzyskanych wyników uwzględniamy także czynnik związany ze złożonością zapro-
jektowanego regulatora.
4. Zaprogramować i zrealizować program badań symulacyjnych na temat roli względnej wagi
czynnika związanego z przebiegiem wielkości sterującej we wskaz
niku jakości regulacji.
Inne parametry niezbędne do syntezy regulatora ustalamy przy tym na poziomie wynikają-
cym z poprzednich badań.
W ramach zadań nadobowiązkowych:
1. Należy przeprowadzić rozszerzony program badań na temat doboru horyzontów predykcji
P i sterowania M. W tym celu przy ustalonej postaci wskaz
nika jakości wielokrotnie
powtarzamy program badań z trzeciego punktu ćwiczenia dla różnych wartości stałej cza-
sowej w modelu obiektu sterowania stanowiącego połączenie członu inercyjnego
z opóz
nieniem. Czy otrzymane rezultaty pozwalają na wypracowanie wskazówek do pro-
jektowania typu: zalecany wybór horyzontu predykcji przy syntezie regulatora na tle rela-
cji do największej stałej czasowej sterowanego procesu?
2. Zaprojektować regulator predykcyjny MPC dla zadanej transmitancji obiektu sterowanego
dobierając parametry syntezy na podstawie wyników przeprowadzonych wcześniej badań,
a następnie porównać rezultaty kompensacji deterministycznego oddziaływania zakłóca-
jącego w klasycznym układzie regulacji i w ramach sterowania predykcyjnego.
SPRAWOZDANIE POWINNO ZAWIERAĆ:
Y
Y Teksty najbardziej istotnych składników oprogramowania stosowanego w trakcie badań.
Y
Y Najbardziej istotne wyniki przeprowadzonych badań symulacyjnych (wskazana jest po-
stać tabelaryczna i graficzna).
21
Y
Y Odpowiedzi na pytania wymienione w programie ćwiczenia.
22
Zestawienie A.
Rekurencyjna estymacja parametrów modelu procesu.
Rekurencyjne sieci działań są w istocie algorytmami działającymi w trybie nadążnym
tj. uaktualniającymi rezultaty obliczeń stopniowo, po każdym nowo-otrzymanym komplecie
danych. Najczęściej chodzi tu o przetwarzanie danych pomiarowych w czasie rzeczywistym.
Uaktualnienie informacji o wyniku estymacji odbywa się wówczas na bieżąco, wraz z nowo
uzyskanym wektorem pomiarów wejść i wyjść badanego obiektu, poprzez wprowadzenie
poprawki do statystycznych ocen szacowanych parametrów:
cest(k) = cest(k-1)+"cest(k)
Wprowadzana poprawka "cest(k) jest zwykle proporcjonalna do różnicy między ak-
tualnym wynikiem pomiaru wielkości wyjściowej y(k) i wynikiem estymacji (predykcji)
wielkości wyjściowej wynikającej z aktualnej w danej fazie obliczeń, oceny parametrów es-
tymowanego modelu: y(k)=c �(u):
yest(k) = cest(k-1)�(u(k));
W rezultacie otrzymujemy wyrażenie ujmujące zasadę obliczeń w postaci ogólnej:
cest(k) = cest(k-1) + K(k)(yest(k)-cest(k-1)�(u(k))).
Poszczególne metody rekurencyjne estymacji parametrów różnią się przede wszystkim spo-
sobem przekształcania błędu predykcji w poprawkę ocen parametrów modelu, tj. sposobem
wyboru wektora K(k). Wektor korekcyjny K(k) może być interpretowany jako wektor wska-
zujący poziom wpływu na poszczególne składowe cest(k) realizacji różnicy wyjścia z aktu-
alnego modelu badanego procesu i jego pomiaru.
Rekurencyjna sieć działań MNK podaje zasady przetwarzania na bieżąco danych po-
miarowych w celu estymacji parametrów procesu dynamicznego w modelu ARX, zwanym
także modelem najmniejszych kwadratów (ang. LS - least squares). Wektor estymowanych
parametrów w modelu ARX składa się ze współczynników wielomianów A(z-1) , B(z-1) .
Przypomnijmy (oznaczenia jak w zestawieniu A / cz.2):
23
-
= - - + � ;
-
Dysponując aktualnymi, w chwili kT, ocenami parametrów tych wielomianów:
cest(kT) = [a1(kT),a2(kT),...,anA(kT),b0(kT),b1(kT),...,bnB(kT)]T;
można uzyskać predykcję wyjścia obiektu w następnym, (k+1)T takcie poprzez wyliczenie:
= uw(k)cest(k).
W wyrażeniu tym przez uu(k) oznaczono wektor:
uw(kT) = [-y(kT),-y((k-1)T),...,-y((k-nA)T),b0((k-d)T),
b1((k-d-1)T),...,bnB((k-nd-nB-1)T)];
Pomiar wyjścia w chwili kT daje z reguły inny wynik (oznaczymy go przez y(k)) niż
jego jednotaktowa predykcja Różnica stanowi podstawę korekcji ocen pa-
rametrów w algorytmie określany zwykle jako RLS (ang.: Recursive Least Squares):
(k)= (k-1)+K(k)(y(k)-uw(k) (k-1)); (0) = co ;
est est est est
Wektor korekcyjny K(k)jest tu określany przedtem za pośrednictwem wyrażenia:
K(k) = P(k)uwT(k) = P(k-1) uwT(k)/ (1 + uw(k)P(k-1)uwT(k));
w którym macierz P(k-1)znana jest już z poprzedniego taktu.
Trzecia składowa algorytmu jest związana z uaktualnieniem macierzy P:
P(k) = [I - K(k) uw(k)] P(k-1);
Na starcie algorytmu przyjmujemy z reguły: P(0)= wI, gdzie: w >> 0.
Warunki zbieżności algorytmu wiążą się z poprawnością przyjętych założeń odnośnie
struktury modelu, z jego stacjonarnością, własnościami sygnałów wejściowego u(k)
i zakłócającego v(k), a także braku korelacji między nimi. Są one analogiczne do warun-
ków poprawności rozwiązania problemu estymacji parametrów w modelu ARX metodą
najmniejszych kwadratów w trybie off line.
Zastosowanie RLS jest skuteczne w sytuacjach ze stosunkowo nieznacznym poziomem
stosunku mocy zakłóceń do sygnału użytecznego w doświadczeniu. Charakteryzują się
24
one w tej sytuacji stosunkowo szybką zbieżnością. W praktyce często korzysta się z wersji
algorytmu uwzględniającej mechanizm zapominania (tj. preferującej informację najbar-
dziej aktualną) w celu uzyskania mechanizmu śledzenia za wolnymi zmianami w czasie
parametrów obserwowanego procesu. Rekurencyjny algorytm metody ważonych NK róż-
ni się nieznacznie sposobem obliczania wektora K(k). Modyfikację tego typu przedsta-
wiono poniżej.
Problem śledzenia za wolnymi zmianami właściwości obserwowanego obiektu może
być rozwiązany przy pomocy jednej z wersji metody ważonych najmniejszych kwadra-
tów, preferującej aktualne wyniki pomiaru kosztem starszych, częściowo zdezaktualizo-
wanych. Efekt tego typu najczęściej osiąga się wprowadzając do bazowego wskaz
nika
kwadratowego  współczynnik zapominania , otrzymując wskaz
nik jakości estymacji
w postaci:
=  - - <  < .
"
=
Waga kwadratu błędu maleje wówczas wykładniczo wraz ze wzrostem liczby taktów
(k i)dzielących dany, i-ty pomiar od momentu aktualnego k. Można wykazać [4], że
dane odległe o więcej niż 3/(1-) taktów od chwili bieżącej posiadają wagę mniejszą
niż5%wagi danych aktualnych.
Uwzględnienie współczynnika zapominania jedynie nieznacznie modyfikuje treść obli-
czeń; w wyrażeniach algorytmu RLS pojawia się współczynnik .Są one przedstawione
nieco dalej.
Pakiet Identyfikacji zawiera wiele m-funkcji przeznaczonych do reku-
rencyjnych metod identyfikacji obiektów dynamicznych różniących się typem stosowanej
struktury modelu. Podstawowe informacje o tym pakiecie oprogramowania zgrupowano
w ramach zestawienia D. W tym miejscu podkreślmy jedynie, że do stosowania modelu
ARX przystosowana jest m-funkcja rarx. Poznając różne wersje realizacji obliczeń przy
jej pomocy możemy m.in. uzupełnić wiedzę nt. metod zapominania informacji przesta-
rzałej.
Rekurencyjna sieć działań rozszerzonej metody NK (ang. RELS: Recursive Ex-
tended Least Squares Method) wiąże się z modelem transmitancyjnym obiektu w strukturze
ARMAX:
25
- -
= - - + � ,
-
równoważnym pseudo-liniowej postaci równania regresyjnego:
y(k)= uw(k) - + �(k),
z wektorem estymowanych parametrów:
c = [a1 , a2 , ... , anA , b0 , b1 , ... , bnB , c1 , ..., cnC]T .
Sekwencja obliczeń jednego taktu rozszerzonej metody NK nie różni się formalnie (meryto-
rycznie, oczywiście, wektory c i uw(k) mają szersze znaczenie ) od prezentowanych wy-
żej, w algorytmie RLS. Przytoczymy je tutaj w wersji ze współczynnikiem zapominania :
cest((k)= cest(k-1)+K(k)(y(k)-U[k]cest(k-1)); cest(0)=co ;
K(k) = P(k-1) uwT (k) / [ + uw(k) P(k-1) uwT(k)];
P(k)=[P(k-1)-P(k-1)uwT(k)uw(k)P(k-1)/[+uw(k)P(k-1)uwT(k)]]/ ;
Wektor zmiennych regresyjnych:
uw(k) = [-y(k-1)- ... -y(k-nA),u(k-d-1),...,u(k-d-nB),
�(k-1) ,..., �(k-nC)];
jest tu na bieżąco estymowany za pośrednictwem różnicy wyjścia z pomiaru i z mo-
delu:
�
Wielomian C(z-1) w modelu ARMAX powinien spełniać warunek pozytywnej realności;
tj. po podstawieniu z = ej�T , wyrażenie:
re{1/C(z-1) -1/2} > 0; dla wszystkich: 0 < �T < Ą.
Inne warunki zbieżności algorytmów RELS są analogiczne, jak RLS.
RELS pracują poprawnie nawet przy dość znacznej roli zakłóceń w modelu procesu. Cha-
rakteryzuje je jednak słabsza zbieżność (niż RLS) w fazie początkowej. Oszacowania
współczynników ci w wielomianie charakteryzującym tor zakłóceń ogólnie rzecz biorąc
zbiegają się wolniej niż oszacowania ai i bi .
26
Rekurencyjna metoda zmiennych instrumentalnych (ang. RIV: Recursive In-
strumental Variables Method) formułuje również sieć działań zbliżoną do RLS:
cest(k)=cest(k-1)+ K(k)(y(k)-uw[k]cest(k-1)); cest(0)=co;
K(k) = P(k-1) wwT (k) / [1 + uw(k) P(k-1) wwT (k)];
P(k)={P(k-1)-P(k-1)wwT(k)uw(k)P(k-1)/[1+uw(k)P(k-1)wwT(k)]}/.
Główna różnica w stosunku do algorytmów przedstawionych wcześniej polega na zastosowa-
niu dwu wektorów zmiennych regresyjnych:
uw(k) = [-y(k-1) - ... -y(k-nA), u(k-d-1),..., u(k-d-nB)];
ww(k) = [-v(k-1) - ... -v(k-nA), u(k-d-1),..., u(k-d-nB)],
z których jeden jest związany z, określanymi na podstawie aktualnego modelu, zmiennymi
instrumentalnymi:
v(k)= uw(k)cest(k-1).
Wektor estymowanych parametrów pozostaje analogiczny jak w modelu ARX. Działania
zgodne z metodą RIV sprzężone są z reguły ze stosowaniem w fazie początkowej, przy
znacznej rozbieżności ocen parametrów i rzeczywistości, klasycznego algorytmu naj-
mniejszych kwadratów RLS.
Stosując algorytmy RIV otrzymujemy z reguły, kosztem nieco większych trudności obli-
czeniowych, istotnie lepszą jakość oszacowań parametrów transmitancji badanego obiek-
tu. Wynika to z niewrażliwości metody zmiennej instrumentalnej na sprzężony
z transmitancją obiektu model filtru w torze zakłóceń.
Inne algorytmy rekurencyjne przetwarzania na bieżąco, w trybie nadążnym danych
pomiarowych w celu estymacji parametrów modelu procesu dynamicznego są również stoso-
wane. Wynika to m.in. z braku uniwersalnego narzędzia obliczeniowego tego typu pasującego
do każdej, potencjalnie napotykanej w praktyce, sytuacji. Wymienimy dwa przykłady.
27
Przy stacjonarności parametrów modelu i wyraz
nych preferencjach dla prostoty obliczeń
użytecznym możne być algorytm aproksymacji stochastycznej (STA). Dla modelu typu
przyjmuje on postać:
= - + - - ;
+ - �
w której zaleca się stosowanie: c1 e" 1, c2 e" 0, c3 > 0, 0 d" � d" 0,5. Zawarto w nim szcze-
gólnie prosty sposób obliczania, malejącego zresztą w czasie do wartości zerowej, wektora
wzmocnienia korekcyjnego. Wadą tego algorytmu jest jego stosunkowo wolna zbieżność.
Zbliżony typ obliczeń, z istotnie innym sposobem wyznaczania wektora wzmocnienia repre-
zentuje algorytm stochastycznego gradientu:
= - + - -
+
Sposób wyliczania wektora korekcyjnego nie jest związany tu z numerem taktu obliczeń.
Np. przedziały dla dobieranych współczynników algorytmu : 0 < c1 < 10;
0 < 0,9 < c2 < 1 <2.
28
Zestawienie B.
Cyfrowe układy regulacji adaptacyjnej.
Cele regulacji adaptacyjnej są analogiczne do celów regulacji klasycznej. W tym przy-
padku chodzi również o zapewnienie układowi sterowania zadanych własności dynamicznych
i własności stanu ustalonego. Występuje tu jednak element dodatkowy, znacznie utrudniający
projektowanie. Elementem tym jest wymaganie utrzymania wymaganych własności układu
regulacji nawet przy nadmiernych (z klasycznego punktu widzenia), lecz stosunkowo wol-
nych zmianach właściwości obiektu sterowania i charakterystyk oddziaływujących nań zakłó-
ceń. W praktyce problem adaptacji sprowadza się do konieczności wprowadzenia w skład
reguły sterowania mechanizmu dopasowania parametrów regulatora do aktualnych własności
sterowanego procesu. Zastosowanie techniki cyfrowej osłabia istotnie ograniczenia natury
sprzętowej pozwalając projektantowi skupić się na aspekcie programowym, tj. problemie
konstrukcji zasad modyfikacji zmian parametrów w stosowanych algorytmach sterowania.
Główna linia podziału przy merytorycznej klasyfikacji układów regulacji adaptacyjnej prze-
biega na tle posiadanych zasobów informacji wyjściowej i sposobu przetwarzania on line in-
formacji bieżącej uzyskiwanej w trakcie ich działania.
Układ z programowalnymi zmianami parametrów regulatora jest stosunkowo
najprostszą wersją układu adaptacyjnego. Z założenia przed jego syntezą wymagane jest jed-
nak posiadanie pełnej informacji na temat zależności właściwości sterowanego obiektu od
określonych, podlegających (bezpośrednim lub pośrednim) pomiarom, tzw. zmiennych wio-
dących. Zmiennymi wiodącymi wybierane są różnorodne wielkości fizyczne. Mogą one być
reprezentowane zarówno przez oddziaływania zewnętrzne (np. prędkość lub wysokość lotu
w przypadku aparatu latającego), jak i przez długość przedziału czasowego który upłynął od
momentu rozpoczęcia procesu technologicznego o charakterze wsadowym, czy też przez po-
łożenie zmieniającego się punktu pracy na charakterystyce nieliniowej sterowanego obiektu.
W fazie projektowania ustala się w tym przypadku zależność funkcyjną (o charakterze deter-
ministycznym lub rozmytym) wiążącą wartości modyfikowanych w czasie rzeczywistym pa-
rametrów regulatora z bieżącą wartością zmiennych wiodących. Podkreślmy, że obwód
adaptacji w układzie z programowalnymi zmianami parametrów regulatora jest w istocie ob-
29
wodem otwartym. Zalety i wady otwartego obwodu sterowania są dobrze znane. Warunki
istnienia zmiennej wiodącej i posiadania pełnej informacji wyjściowej ograniczają poważnie
liczbę zastosowań.
W układach pośredniej regulacji adaptacyjnej konieczne jest posiadanie informa-
cji wyjściowej na temat struktury modelu sterowanego obiektu, umożliwiającej prowadzenie
procesu identyfikacji w trybie nadążnym. Uzyskiwane w wyniku obliczeń rekurencyjnych
oceny parametrów tego modelu stanowią dane wyjściowe do syntezy, także w trybie on line,
parametrów regulatora.
Utrzymując się w ramach tej ogólnej reguły możemy stosować, dopasowane do modelu
obiektu i zakłóceń w rozwiązywanym zadaniu sterowania, różnorodne rekurencyjne algoryt-
my identyfikacji i metody syntezy regulatora. Stosunkowo często podstawę rozwiązania pro-
blemu adaptacji stanowią np. znane z poprzednich ćwiczeń:
algorytmy estymacji parametrów RLS i RELS dostosowane do modeli dynamiki sterowa-
nego obiektu (odpowiednio typu ARX i CARMA),
zasady regulacji minimalnowariancyjnej (regulator MVC w przypadku silnych oddziały-
wań zakłócających), czy też ustalenie położenia biegunów, lub transmitancji wzorcowej
układu wynikowego (regulator PPC, lub PZPC w sytuacjach związanych ze stosunkowo
niskim poziomem zakłóceń) przy syntezie regulatora.
Strukturę rozwiązania zadania adaptacji pośredniej podajemy przekazując obie główne
składowe rozwiązania, np. RLS / MVC, RELS / PPC lub RLS / PZPC. Nie wyklucza to,
oczywiście, stosowania innych, niesygnalizowanych tu metod. Można wymienić także np.
połączenia sieci działań w postaci STO/PID, RIV/SC, na podstawie bieżących wyników
identyfikacji stosuje się reguły regulacji predykcyjnej itp. Czynnikiem decydującym
o wyborze struktury algorytmicznej adaptacji jest osiągnięcie zadawalających wskaz
ników
jakości rozwiązania danego zadania sterowania. Podkreślmy, że oba główne elementy struk-
tury algorytmicznej są niezależne; sposób i kryteria syntezy regulatora nie są związane
z przyjętą metodą identyfikacji. W trakcie syntezy regulatora, na miejsce parametrów modelu
sterowanego obiektu podstawiane są ich aktualne oceny.
30
Podkreślmy także, że zadaniem regulatora w układzie regulacji adaptacyjnej jest wygene-
rowanie sygnału sterującego, który oprócz swej klasycznej roli, ma ułatwić proces bieżącej
identyfikacji obiektu sterowania. Problem ten jest istotny m. in. z tego względu, że proces ten
odbywa się z natury rzeczy w zamkniętej strukturze układu regulacji, implikującej sprzężenia
korelacyjne ciągów czasowych oddziaływań sterujących i wielkości regulowanych. Najprost-
szy sposób zagwarantowania identyfikowalności modelu procesu sterowania związany jest
z zastosowaniem zewnętrznego sygnału trwale pobudzającego dostatecznego rzędu i dopusz-
czalnej amplitudy.
Układy bezpośredniej regulacji adaptacyjnej budowane są z pominięciem fazy
identyfikacji modelu sterowanego procesu. Wdrażany jest tu algorytm rekurencyjny najmniej-
szych kwadratów realizujący zasadę bezpośredniej identyfikacji współczynników wielomia-
nów R, S, T z równania regulatora [4] na podstawie reprezentowanych w nim realizacji
ciągów czasowych pomiarów wielkości regulowanej y(k), znajomości sygnału odniesienia
w(k), generowanych sterowań u(k)i znajomości przebiegów uogólnionej, skonstruowanej
z uwzględnieniem zadania sterowania wielkości wyjściowej Y*(k).
Istota bezpośredniej regulacji adaptacyjnej tkwi tu więc w konstrukcji umożliwiającego
ten proces regresyjnego modelu predykcyjnego typu:
Y*(k) = Rposr u(k-d) + Sposr y(k-d)  Tposr w(k-d) + �(k).
W modelu tym:
Zmienne regresyjne stanowią aktualne (z uwzględnieniem opóz
nienia d reprezentowa-
nego w modelu obiektu) składowe wektorów ciągów czasowych y(k), u(k), w(k)
z modelu regulatora.
Wektor estymowanych parametrów powstaje w wyniku wspólnego ujęcia współczynni-
ków wielomianów pośredniczących Rposr, Sposr, Tposr związanych prostą (najczęściej
poprzez wspólny składnik multiplikatywny) zależnością z wielomianami: R, S, T z
równania regulatora.
Przy pomocy Y*(k) oznaczono tu tzw. uogólnioną wielkość wyjściową modelu predyk-
cyjnego zależną zarówno od wielkości interpretowanych tu jako zmienne regresyjne, jak
i od stosowanego prawa regulacji. Sposób wyrażenia tego składnika jest sprawą kluczową
w trakcie rozwiązania problemu sterowania.
31
 Czynnik �(k)ujmuje w równaniu oddziaływujące zakłócenia.
Przytoczymy, za [4], kilka konkretnych przykładów modeli tego typu dla modelu stero-
wanego obiektu w postaci: A y(k) = B u(k-d).
W celu samo nastrajania regulatora na zadane położenie biegunów układu regulacji (rów-
noważne z ustaleniem wielomianu charakterystycznego A*(z-1)) powstaje model predyk-
cyjny w postaci:
A* Apom y(k) = R B u(k-d) + S B y(k-d) + �(k).
Rolę wielomianu pomocniczego Apom omówiono ( PP- regulacja ) w jednym z zestawień
z rozdz.2. Obecność wielomianu B(z-1) istotnie komplikuje proces dostrajania parame-
trów regulatora zwiększając rząd stosowanego modelu i wprowadzając dodatkowe dzia-
łania natury obliczeniowej.
Przy regulacji adaptacyjnej z ustalaniem położenia biegunów i zer (PZP1) uzyskujemy
natomiast model znacznie prostszy:
A* Apom y(k) = R u(k-d) + S y(k-d) + �(k).
Przy celowości zastosowania zasad regulacji minimalnowariancyjnej wprowadza się mo-
del obiektu z uwzględnieniem toru zakłóceń probabilistycznych (np. poprzez dodatkowy
wielomian C(z-1) w modelu ARMAX). Komplikuje to, oczywiście, i model predykcyj-
ny, który przyjmuje nieco bardziej złożoną postać:
y(k) = R u(k-d) + S y(k-d) + F �(k) + (C  1)[ F �(k)-y(k)].
Po zdefiniowaniu równania regresyjnego konstruowana jest rekurencyjna sieć działań ce-
lem której jest początkowe dostrojenie i następnie śledzenie w trakcie pracy ewentualnych
zmian położenia zestawu parametrów regulatora R, S, T optymalnych z punktu widzenia
wskaz
nika kwadratowego:
J = min Łkj=1 ąk-j[(Y*(k)-Rposru(k-d)-Sposry(k-d)+Tposrw(k-d))]2.
32
Zdolność nadążania algorytmu za ewentualnymi zmianami optymalnych nastaw zapewnia
zastosowanie techniki współczynnika zapominania ą<1.
Podkreślmy, że rozpoczęcie pracy algorytmu wiąże się z wprowadzeniem warunków po-
czątkowych typu: {R, S, T]=0; P(0)= �I, gdzie �>>0, w celu osiągnięcia okre-
śloności rozwiązania zmienne regresyjne powinny mieć charakter  trwale pobudzający .
Nadrzędna warstwa diagnostyki i nadzoru: odgrywa podstawową rolę w systemie
automatyki. Zadaniem umiejscawianych tu działań jest m.in. ocena jakości istniejących
w ramach systemu mechanizmów adaptacji. Kontroli podlega zarówno proces identyfikacji
jak i jakość nastaw regulatora. Mogą tu być umieszczone algorytmy umożliwiające udzielenie
odpowiedzi na pytania typu:
" Czy wyjściowo wybrana struktura modelu sterowanego procesu jest właściwa?
" Czy uwzględniane w trakcie estymacji pobudzenia są dostatecznie urozmaicone?
" Czy zastosowana metoda zapewnienia zdolności śledzenia zmian parametrów modeli (np.
poprzez ustalenie współczynnika zapominania) jest racjonalną w świetle napływających
danych bieżących?
" Czy aktualnie obserwowane (krótko  i długoterminowe) tendencje zmian wskaz
ników
oceny jakości przebiegu sterowanych procesów są uzasadnione zmianami czynników ze-
wnętrznych?
Przypomnijmy, że w klasycznej, warstwowej strukturze systemów sterowania możemy
wyróżnić kilka kolejnych warstw sprzętowo  programowych.
1. W ramach warstwy instrumentalizacji i regulacji bezpośredniej umieszcza się m.in.
podstawowe algorytmy przetwarzania sygnałów, jak filtracja, korekcja, algorytmy regu-
latorów, itp.
2. Do warstwy regulacji nadrzędnej zalicza się działania związane z wygenerowaniem
trajektorii odniesienia, kontrolę ograniczeń z alarmowaniem, algorytmy realizacji zmian
punktów pracy, itp.
3. W zadaniach rozwiązywanych w warstwie optymalizacji sterowania i w warstwie adap-
tacji chodzi o zabezpieczenie nadążnego przetwarzanie zmian systemowych (tj. zmian
właściwości obserwowanych oddziaływań i postaci stosowanych modeli procesów) na
33
odpowiednie, oceniane przez pryzmat określonych wskaz
ników jakości, zmiany parame-
trów, a nawet czasem struktur programowych w warstwach sterowania.
4. Do warstwy nadzorującej systemu sterowania (ang.  supervisory control layer ) zali-
cza się bardzo szeroką listę zadań z których oczywiście nie wszystkie muszą być realizo-
wane. Do ważniejszych z nich należą niewątpliwie algorytmy wielowymiarowej analizy
statystycznej sterowanych procesów (ang. MSPC  multivariate statistical process con-
trol) pozwalające m.in. na ustalanie w trybie nadążnym statystycznych niezgodności ist-
niejących reżimów technologicznych w porównaniu z zadanymi w celu oceny na jakim
poziomie jakościowym przebiega aktualnie proces technologiczny. W rezultacie uzyskuje
się monitorowanie jakości połączone z pozyskiwaniem materiału do bazy wiedzy nt. ste-
rowanego procesu. Inną bardzo istotną grupę zadań w tej warstwie realizują także algo-
rytmy detekcji, zlokalizowania i izolacji usuwania uszkodzeń (ang. FDI  faults detection
and isolation) dając możliwość projektowania systemów sterowania odpornych na (nie-
które) uszkodzenia (ang. FTC  fault tolerant control).
34
Zestawienie C.
Sterowanie predykcyjne.
Zasada sterowania predykcyjnego wynika bezpośrednio z celu stosowania regulacji
automatycznej.
W momencie k należy oddziaływać na sterowany obiekt takim sterowaniem u(k), by
predykcja wielkości regulowanej na takty następne: yP(k+j|k); j=1,2,...,P poczy-
niona przy uwzględnieniu stanu obiektu w momencie k - tym przyjęła znane wartości pożą-
dane: y*(k+j|k).
Załóżmy, że obiekt sterowany opisywany jest znanym równaniem różnicowym:
y(k)=-a1y(k-1)-a2y(k-2)-...-any(k-n)+b1u(k-1)+...+bmu(k-m)+�(k);
w którym czynnik �(k) uwzględnia nieznane oddziaływanie zakłócające o zerowej wartości
średniej. W sytuacji, gdy równanie to jest znane, możemy posłużyć się nim by wyliczyć ra-
cjonalną predykcję wartości wielkości wyjściowej: yP(k+1|k) na takt następny przy in-
formacji dostępnej w chwili bieżącej.
yP(k+1|k)=-a1y(k)-a2y(k-1)-...-any(k+1-n)+b1u(k)+...+bmu(k+1-m);
Istotę sprawy stanowi ten fakt, że z posiadania informacji o modelu obiektu wynika rów-
nież możliwość ustalenia wartości jego wielkości wyjściowej na pożądanym (tj.
y*(k+1|k))poziomie. Wystarczy w tym celu jedynie wyznaczyć i zastosować wynikające
z modelu oddziaływanie sterujące:
u(k)=[y*(k+1|k)+a1y(k)+..+any(k-n+1)+b2u(k-1)+..+bmu(k+1-m)]/b1.
W praktyce należy zadbać także o to, by oddziaływanie to nie okazało się zbyt duże,
przekraczające możliwości urządzenia wykonawczego w układzie regulacji. Podkreślmy tak-
że, że samo stosowanie liniowego przybliżenia modelu sterowanego obiektu, uwarunkowane
jest często stosowaniem umiarkowanych przyrostów oddziaływań sterujących względem
punktu pracy.
35
Istnienie ograniczeń i dynamika sterowanego procesu powoduje celowość przeprowadze-
nia dodatkowych zabiegów. Najbardziej istotny postulat dodatkowy zakłada przyjęcie stop-
niowego, wielotaktowego przechodzenia do pożądanej wartości zmiennej regulowanej (np.
zadany poziom w układzie stabilizacji). Wystąpi wówczas konieczność rozwiązywania (rów-
nolegle z zadaniem wyliczania u(k)) sprzężonego z tym zadania określania,  przyszłościo-
wej sekwencji wyjść układu długości P:
Sekwencja tego typu nazywana jest trajektorią odniesienia, lub pożądaną trajektorią ru-
chu sterowanego obiektu, a występującą w nim liczbę naturalną P określa się jako horyzont
predykcji. Uwzględniana przy obliczeniach trajektorii odniesienia liczba M aktywnych (tj.
niezależnie w każdym takcie obliczeń kształtowanych) sterowań nazywana jest z kolei hory-
zontem sterowania.
Trajektorię odniesienia uzyskuje się stosunkowo często na podstawie przyjęcia określo-
nej, możliwie nieskomplikowanej postaci transmitancji (wzorcowej): H*(z-1) układu wyni-
kowego. W roli przykładów można tu wymienić:
przyjęcia za pożądaną trajektorię zachowanie układu o transmitancji równoważnej
z unormowaną transmitancją (przez K0 oznaczono wzmocnienie) sterowanego procesu:
H*(z-1)= G0(z-1)/K0;
zastosowanie w roli generatora wzorcowej trajektorii członu inercyjnego z opóz
nieniem:
H*(z-1)= (1-�) z-1/(1-�z-1 );
zatrzymanie się na prostej wersji czystego opóz
nienia: H*(z-1)= z-r .
Pożądane zachowanie się projektowanego układu wynika wówczas z cech przyjmowane-
go modelu generującego w trybie nadążnym trajektorie odniesienia. Powstaje wówczas wy-
raz
ne podobieństwo metodyki sterowania predykcyjnego do syntezy regulatora cyfrowego
z ustaleniem zer i biegunów wynikowego układu regulacji. Podkreślmy, że koncepcja trajek-
torii odniesienia jest jednak ogólniejszą od koncepcji modelu odniesienia określającego wzor-
cowe własności dynamiczne projektowanego układu..
36
Innym zabiegiem osłabiającym wymagania odnośnie poszczególnych elementów ciągu
sterowań jest zastosowanie w trakcie ich wyznaczania optymalizacji zorientowanej na mini-
malizację wskaz
nika kwadratowego jakości regulacji opartego najczęściej na sumie ważonych
kwadratów elementów ciągu sterowań (lub, zamiennie, jednotaktowych przyrostów tych ele-
mentów) i kwadratów odchyłek uzyskiwanej trajektorii od pożądanej.
Zadanie syntezy sterowania predykcyjnego polega na znalezieniu takiego ciągu
sterowań {up(k|k), up(k+1|k),... up(k+P-1|k)}, który minimalizuje skonkrety-
zowaną wersję wskaz
nika jakości regulacji J w ogólnej postaci przedstawianego wyrażeniem:
J = Ł [(yP(k+j|k)- y*(k+j|k))2 + � �u2(k+j|k)};
w którym sumowanie odbywa się po wszystkich numerach taktów horyzontu predykcji: j =
1,...,P, a symbol � sygnalizuje możliwość uprzedniego przeprowadzenia określonej
operacji (np. odjęcie stałej, różnicowanie) na elementach ciągu czasowego sterowań.
Przyjmuje się przy tym następujące założenia:
Znany jest model obiektu umożliwiający predykcję P kolejnych wartości wyjścia obiektu:
yp(k) ={yp(k+1|k),yp(k+2|k),... yp(k+P|k)},
tj. aż do (k+P)-tego taktu włącznie.
Znane są na bieżąco, niezbędne do tej predykcji sekwencje poprzednich (tj. poprzedzają-
cych moment aktualny) sterowań:
{u(k-1),u(k-2),... u(k-m~+1)}
i (ewentualnie) wielkości wyjściowej:
{y(k),y(k-1),... y(k-n~+1)}
sterowanego obiektu.
Dla przedziału czasowego [k, k+P] generowana jest, w wyniku zastosowania niezależ-
nie sformułowanej reguły, pożądana trajektoria ruchu obiektu sterowania.
Istnieje możliwość skonkretyzowania tego sformułowania w wielu wersjach, różniących
się potencjalnie:
37
 przyjętą postacią wskaz
nika jakości regulacji  znane są tu rozwiązania począwszy od
minimalnowariancyjnych, w ramach których we wskaz
niku wyrażana jest wyłącznie skła-
dowa związana z dokładnością odtworzenia trajektorii, aż po wskaz
niki związane z mini-
malizacją energii w których istotna jest wyłącznie energia sterowania  najczęściej jednak
przyjmuje się rozwiązania w których występują oba składniki;
stosowanym modelem sterowanego obiektu dynamicznego  główna linia podziału prze-
biega od modeli nieparametrycznych opartych na charakterystykach czasowych obiektu
do modeli parametrycznych, tworzonych na podstawie równań różnicowych; w niektó-
rych wersjach modelu sterowanego procesu uwzględnia się również model toru zakłóceń;
sposobem wypracowania trajektorii odniesienia  typowym chwytem jest wspomniane już
wyżej automatyczne generowanie trajektorii odniesienia w postaci równoważnej poten-
cjalnej trajektorii obiektu o ustalonej transmitancji.
W związku z tym możliwe jest sformułowanie i rozwiązanie w postaci ogólnej różnych
wersji zadania syntezy sterowania predykcyjnego. Dodajmy, że rozmaitość badanych i stoso-
wanych w praktyce rozwiązań wzrasta jeszcze przy potencjalnej konieczności uwzględnienia
ścisłych ograniczeń dopuszczalnej amplitudy oddziaływań sterujących. Zadania tego typu
ujmuje teoria układów predykcyjnych z ograniczeniami dopuszczalnych sterowań lub teoria
nieliniowych układów predykcyjnych.
Nieznajomość modelu sterowanego obiektu prowadzi do układów typu adaptacyjnego
(APC  ang. Adaptive Predictive Control). Należy przeprowadzać bowiem wówczas proces
identyfikacji w trybie nadążnym i posługiwać się, uzyskanym w trakcie tego procesu, mode-
lem opartym na ocenach liczb określających jego strukturę: n, m ( tj. stopnie wielomianów
licznika i mianownika transmitancji - rząd modelu), wielkość opóz
nienia di estymaty wystę-
pujących w nim współczynników: ai, bj.
Charakterystyczną wersją regulacji predykcyjnej jest badana w ćwiczeniu wersja z mo-
delem odpowiedzi impulsowej nazywana jest MPC (ang. model predictive control). Charakte-
rystyka impulsowa sterowanego procesu jest, jak wiemy, stosunkowo łatwa do uzyskania
doświadczalnie. W wersji tej do obliczeń trajektorii stosowany jest model obiektu wynikający
z dyskretnego równania splotu. Przypomnijmy, że wyrażenie splotu w wersji dyskretnej:
38
y(k)= g(0)u(k)+g(1)u(k-1)+...+g(i)u(k-i)+...+ �(k);
po ograniczeniu ważącej wpływ sterowań charakterystyki impulsowej g(k)do liczby mwy-
razów, staje się z
ródłem przybliżonego modelu liniowego obiektu dynamicznego. Ucięcie
liczby wyrazów g(k) jest równoważne z przyjęciem założenia o nieistotności wyrazów
póz
niejszych tej charakterystyki, nazywanej często także funkcją wagową sterowanego
obiektu: |g(i)| H" 0 dla i > m. Przy racjonalnym wyborze m model ten może umoż-
liwiać, podobnie jak model w postaci transmitancji, dostatecznie dobrą predykcję wielkości
wyjściowej.
W ramach założeń wyjściowych regulacji MPC przyjmuje się także kwadratowy, dwuskład-
nikowy wskaz
nik jakości zależny od wyliczanych predykcji wielkości wyjściowej obiektu,
trajektorii odniesienia i sterowań (lub przyrostów sterowań). W sytuacji najprostszej z obiek-
tem SISO rozwiązanie stanowi wektor wzmocnień KMPC różnic między pożądaną wartością
wektora trajektorii odniesienia w przedziale czasowym równym horyzontowi predykcji:
W = [y*(k+1|k), y*(k+2|k), ..., y*(k+P|k)];
i wektorem trajektorii predykcji wielkości wyjściowej obiektu regulacji:
yp = [yp(k+1|k),yp(k+2|k),... yp(k+P|k)],
uzyskiwanym przy założeniu zerowych przyrostów sygnału sterującego poza strefą określaną
przez horyzont sterowania.
Na zakończenie przedstawimy twierdzenie o rozwiązaniu tak postawionego zadania.
Dowód możemy znalez
ć w literaturze ([4]).
Twierdzenie. Przyrosty wielkości wyjściowej regulatora minimalizującego wartość wskaz
ni-
ka:
J = Ł [(yP(k+j|k)- y*(k+j|k))2 + � "u2(k+j|k)};
w którym sumowanie odbywa się po wszystkich taktach horyzontu predykcji: j=1,...,P;
a w końcowej części: j=M+1,...,P; tego przedziału przyrosty sterowania są zerowe,
można wyznaczyć z wyrażenia:
"u(k)= u(k)  u(k-1) = KMPC ( y*(k) - yp(k)).
39
W wyrażeniu tym y*(k) i yp(k) są wektorami utworzonymi z odpowiednich fragmen-
tów trajektorii odniesienia i trajektorii predykcji wielkości wyjściowej sterowanego
obiektu:
y*(k) = [ y*(k+1|k), y*(k+2|k),..., y*(k+P|k)]T;
yp(k) =[yp(k+1|k), yp(k+2|k),..., yp(k+P|k)]T;
KMPC jest natomiast wektorem współczynników wzmocnienia regulatora predykcyjnego.
KMPC obliczany jest jako pierwszy wiersz wyrażenia macierzowego:
[QTQ + � I ] QT ;
w którym Ijest macierzą jednostkową, a Q powstaje w wyniku rozmieszczenia elementów
funkcji wagowej sterowanego obiektu g(k)w następującym układzie:
�ł łł
�ł śł
�ł śł
�ł śł
�ł śł
=
�ł - - śł
�ł śł
�ł śł
�ł śł
�ł - - - śł
�ł �ł
40
Zestawienie D.
Pakiet oprogramowania narzędziowego .
Pakiet identyfikacji zawiera około siedemdziesięciu m-funkcji pogrupowanych w kilku
zestawach tematycznych. Celem danego zestawienia jest zwięzłe przedstawienie jedynie naj-
bardziej podstawowych informacji, niezbędnych do zapoznania się z możliwościami tego
pakietu w zakresie realizacji części obliczeń w ramach ćwiczeń tego rozdziału.
Tablica o nazwie   jest stosowana w pakiecie  Identyfikacja do przed-
stawienia w postaci standardowej modelu procesu dynamicznego. Ogólna struktura tego mo-
delu mieści się w ramach równania:
A(z-1) y(k) = [B(z-1)/F(z-1)] u(k-nk) + [C(z-1)/D(z-1)] �(k).
Przez A,B,C,D,F oznaczono w tym równaniu wielomiany (stopnia: na,... ,nf) opera-
tora opóz
nienia o 1 takt z-1, nk jest tu wartością opóz
nienia sprzężoną ze sterowanym wej-
ściem u(k), �(k) reprezentuje natomiast szereg czasowy typu białego szumu. Zauważmy,
że w ramach tego równania (przy założeniu: u(i) a" 0) mieszczą się także modele parame-
tryczne ciągów czasowych.
W tablicy  theta można umieścić również model liniowego obiektu dynamicznego
w wersji ciągłej, z wielomianami operatora  s Laplace a, a także model obiektu wielowy-
miarowego. Model w postaci standartowej najprościej można uzyskać przy pomocy
m-funkcji:   po uprzednim zdefiniowaniu wielomianów w równaniu różnicowym.
W przypadku modelu ARMAX stosujemy np. sekwencję podstawień typu:
O szczegółach reprezentacji modelu w formie  theta dowiadujemy się poprzez wywołanie
  . Wywołanie m- funkcji: wyświetla natomiast wszyst-
kie informacje zawarte w tablicy o nazwie   .
41
Tablica  theta zawiera:
W pierwszym wierszu informacje o estymacie wariancji szumu � o czasie próbko-
wania, o liczbie sterowanych wejść ( obiektu i (kolejno) o stopniach stosowanych w
modelu wielomianów:
W wierszu drugim znajduje się informacja o wskaz
niku rozbieżności: dane pomiarowe-
wyjście z modelu (ang. FPE: Final Prediction Error) i dacie wygenerowania modelu.
W wierszu trzecim tablicy umieszczono ( w postaci uproszczonej, bez oczywistych jedy-
nek i zer ) wektory estymowanych współczynników wielomianów A,B,C,D i F.
Następne wiersze tablicy   przekazują (kolejno) wynikową ocenę macierzy kowa-
riancji wszystkich ocen parametrów w modelu.
Proces symulacji doświadczenia możemy oprzeć na modelu typu   i deklara-
cji stosowanych oddziaływań wejściowych: sterujących i zakłócających. Przykładowa se-
kwencja instrukcji:
Dyskretny model theta (o nazwie  mdt ) może być przetransformowany:
mct = thd2thc(mdt); do modelu z czasem ciągłym;
[A,B]= th2arx (mdt); do modelu ARX;
[g, phiv]= th2ff(mdt); do modelu charakterystyk widmowych;
[A,B,C,D,F]= th2poly(mdt); do ocen współczynników wielomianów;
[A,B,C,D,K,X0]= th2ss(mdt) ; do modelu zmiennych stanu;
[licz, mian]= th2tf (mdt) ; do modelu transmitancyjnego;
[zb, k]= th2zp(mdt); do modelu o jawnie wyrażonych zerach i biegunach.
Przypomnijmy, że w  Matlabie istnieją mechanizmy zapewniające możliwość transformacji
pomiędzy każdą parą ogólnie przyjętymi form opisu dynamiki obiektów liniowych.
Odfiltrowanie wysokoczęstotliwościowych składników zakłócających z sygnału  z moż-
na uzyskać wprowadzając filtr dolnoprzepustowy za pośrednictwem instrukcji:
42
zf=ifilt(z,rzfiltru,omega). Argument wejściowy  rzfiltru określa rząd,
natomiast  omega częstotliwość graniczną automatycznie dobieranego filtru Butterwortha.
Usunięcie składowej stałej, lub liniowej tendencji do powolnych zmian można
dokonać przy pomocy: zd=dtrend(z,rztrendu,podzielenie) . Możliwe jest
również, biorąc pod uwagę sens argumentów wejściowych usuwanie tendencji w określo-
nych, określonych przez wektor  podzielenie przedziałach.
Wybór struktury modelu obiektu dynamicznego, przy założeniu jego liniowości, spro-
wadza się głównie do określenia stopnia wielomianów reprezentowanych w transmitancjach
torów sterowania i zakłóceń. W ramach pakietu  Identyfikacji wybór rzędu modelu
sprowadza się do przeliczenia wszystkich potencjalnie możliwych, w świetle posiadanej in-
formacji, wersji rozwiązania, a następnie wyboru najlepszej z nich z punktu widzenia okre-
ślonego wskaz
nika jakości. Odbywa się to w ramach struktury ARX; najprostszej, a przez to
najmniej angażującej moc przeliczeniową komputera przy pomocy następujących pomocni-
czych m-funkcji:
NN= struc(nA,nB,nk) - jest instrukcją pomocniczą; generuje macierz NN  wyli-
czającą wierszami wszystkie, deklarowane przy pomocy argumentów wejściowych, np.:
nA = 2 : 5, nB = 2 : 5, nk = 3:5,
kombinacje standartowego wektora rzędów wielomianów A,B i opóz
nienia d w równa-
niu różnicowym o strukturze ARX: nn = [na, nb, nk].
V = arxstruc(z,zk,NN) - wylicza optymalne parametry i odpowiadające im
wartości kwadratowej funkcji start modeli ARX we wszystkich wersjach nn ( ujętych w
ramach macierzy  NN  ). Istnieje możliwość podziału danych z pomiarów: z = [u,
e] na dwa podzbiory: podstawowy  z i kontrolny  zk .
V = ivstruc(z,zk,NN) - spełnia rolę analogiczną do procedury poprzedniej stosu-
jąc technikę zmiennych instrumentalnych.
nn = selstruc(V,c)- procedura selekcyjna, wyznaczająca na podstawie rezultatów
działania:  ,  optymalny, w sensie wybranego kryterium c,
rząd modelu. Informacje odnośnie możliwości i techniki wyboru stosowanego przy tym
kryterium znajdują się w tekście pomocy tej m-funkcji.
43
Przytoczymy przykład typowej sekwencji instrukcji prowadzącej do wyboru rzędu modelu:
NN = struc(2,2,1:3); % na pocz" tek ustalmy wielko" " opó" nienia;
V=arxstruc(z(1:200,:),z(201:400,:),NN); %dane podzielono na 2 zbiory
nn = selstruc (V,0); nk = nn(3) % opó" nienia zosta" o ustalone;
NN = struc (1:4, 1:4, nk); % teraz ustala si" rz" dy wielomianów A,B;
va = arxstruc (z(1:200,:), z(201:400,:), NN);
vi = ivstruc (z(1:200,:), z(201:400,:), NN);
nna = selstruc (va); nni = selstruc (vi);
Estymacja parametrów modelu przeprowadzana jest w pakiecie  Identyfikacji
bezpośrednio na podstawie nagromadzonych w trakcie pomiarów danych dla następujących
modeli ( w opisie przedstawiono jedynie najprostsze wersje argumentów wejściowych m-
funkcji):
th = ar (y, n) - AR model sygnału stochastycznego n- tego rzędu;
th=armax(z,nn)-stosuje się dla modelu ciągu ARMA lub modelu układu ARMAX;
th = arx (z, nn) - estymator MNK dla modeli AR, lub ARX;
th = bj(n,nn) - model obiektu Boxa-Jenkinsa, określony przez 4 wielomiany;
th = ivar(y,na)- model ARX estymowany metodą zmiennych instrumentalnych;
th = iv4(z,nn) - estymacja z zastosowaniem zmiennych instrumentalnych;
th = oe(z,nn) - estymacja parametrów modelu dynamiki typu  output-error ;
th = pem(z,nn)- model predykcji PEM stosowany do układów wielowejściowych.
Dodajmy, że dla wszystkich wyżej wymienionych modeli parametrycznych potencjalnie
możliwa jest także rekurencyjna estymacja parametrów. Odpowiednie funkcje mają analo-
giczne do wymienionych nazwy poprzedzone literą  r (
, itd.).
Modele nieparametryczne mogą być wyznaczane przy pomocy:
R = cov(z,M)- wyliczenie (M-1)wartości macierzy kowariancyjnej sygnałów, dane
o przebiegach których umieszczono w kolejnych kolumnach argumentu z.
44
 [ir,R,cl] = cra(z,M,na) - estymacja odpowiedzi impulsowej (ir) badanego
obiektu przy pomocy analizy korelacyjnej sygnałów (z = [y u]); wyznacza się tak-
że dla niej szerokość przedziału ufności (cl, przy: ą = 0,01 ) i (dodatkowo) macie-
rze kowariancyjną sygnału R ;
[G,fi] = spa(z,M,w) - estymacja empirycznej transmitancji częstotliwościowej
systemu G(exp(j�)) dla � = ([1:N]/N) " (Ą/T); łącznie z analizą spek-
tralną Ćv(�) zakłóceń w modelu OE;
Przykład sekwencji instrukcji z zastosowaniem analizy nieparametrycznej:
Proces weryfikacji uzyskanego modelu pozwala ocenić, czy uzyskany w wyniku
działań w poprzednich etapach identyfikacji model w dostatecznym stopniu odzwierciedla
właściwości badanego obiektu. Wymienimy ważniejsze funkcje bezpośrednio ułatwiające
ocenę przeprowadzonego procesu identyfikacji:
[e,r] = resid(z,th) - obliczenia poszczególnych residuów e, funkcji autoko-
relacyjnej ciągu residuów i funkcji korelacji wzajemnej sterowania i residuów (e,u);
compare (z,th,k) - przeprowadza proces symulacji z predykcją (na k- taktów
wprzód) wyjścia yh(t) dla uzyskanego modelu przy zadanych oddziaływaniach
u(i); rezultat jest wykreślany łącznie z bezpośrednimi pomiarami wyjścia y;
idsimsd (u,theta,N,szum) - na podstawie zadeklarowanego u i wyników
identyfikacji (model theta) generuje się N, losowo wybranych modeli, ze zbioru ufno-
ści wektora parametrów theta. Następnie wylicza się i odtwarza na monitorze prze-
biegi y(k). Jeżeli szum = `noise' do przebiegów dodawane są zakłócenia.
Procedury prezentacji wyników estymacji stanowią istotne uzupełnienie pakietu
i spełniają istotną rolę przy weryfikacji uzyskanego modelu. Są to w zasadzie nieznacznie
zmodyfikowane funkcje bazowe Matlaba. Wymienimy najważniejsze z nich:
45
 bodeplot ([gl,...gn]) - odtwarza charakterystyki częstotliwościowe w formie
diagramów Bode go. Charakterystyki różnych modeli można porównywać na jednym
wykresie. Argument typu gi (zmienna standartowa ) może być generowany
przez szereg funkcji, np. za pośrednictwem przedstawionych już wcześniej:  ,
 ;
ffplot([gl,...gn]) -charakterystyki częstotliwościowe z liniową skalą;
idplot(z)- wspólne wykresy czasowe wejścia i wyjścia obiektu: z = [y u];
nyqplot([gl,...gn])- diagramy Nyquista dla kilku modeli ;
46
Zestawienie E.
Pakiet oprogramowania sterowania predykcyjnego  mpccmds .
W pakiecie oprogramowania  Model Predictive Control Toolbox , określa-
nym w skrócie poprzez  mpccmds (i sprzężonym z nim  mpcdemos ) zgrupowano m-
funkcje Matlaba przeznaczone do rozwiązywania problemów syntezy, symulacji i analizy
wybranych wersji rozwiązań sterowania predykcyjnego. W zestawieniu przedstawimy niektó-
re z nich, najważniejsze z punktu widzenia zadań realizowanych w trakcie zajęć w kompute-
rowym Laboratorium Metod Sterowania.
Pakiet zawiera około siedemdziesięciu m- funkcji (patrz  help mcpcmds ) zgrupo-
wanych w kilku tematycznie zorientowanych blokach tematycznych.
W pakiecie MPC stosowany jest specyficzny zapis modelu dynamiki w tzw.  MPC step
format powstający przy zastosowaniu uciętej charakterystyki skokowej. Zapis w postaci
 MPC tf format zawiera m.in. informacje na temat transmitancji obiektu. Przekształce-
nia modeli innych typów do tej specyficznej postaci zapisu można przeprowadzić automa-
tycznie korzystając z następujących m- funkcji pakietu spośród ujętych w bloku tematycznym
 model conversion :
poly2tfd  przekształca transmitancję (podaną za pośrednictwem wielomianów w za-
pisie Matlaba) do modelu w postaci  MPC tf format ;
ss2step - przekształca równania stanu (podane za pośrednictwem występujących
w nich macierzy zapisie Matlaba) do modelu w postaci  MPC step format ;
tfd2step - przekształca model przedstawiony w postaci  MPC tf do postaci stoso-
wanej w ramach  MPC step format .
Uwzględniając istnienie bogatego uzupełnienia w postaci zestawu m-funkcji transforma-
cji modeli obiektów dynamicznych w ramach tematycznie sprzężonego pakietu Control
Matlaba (ss2tf, ss2zp, c2d, d2c, itd.) można uznać, że wymienione m- funk-
cje stanowią dostateczne ogniwo konwersji każdego rozpowszechnionego modelu liniowego
dynamiki do postaci wymaganej przez pakiet MPC.
47
Dysponując modelem sterowanego obiektu w postaci MPC step format (argument
wejściowy o nazwie  model ) możemy m. in. zrealizować następujące zadania:
Y
Y Stosując:
 Kmpc = mpccon (model,ywt,uwt,M,P)
możemy obliczyć zestaw współczynników wzmocnień regulatora, stanowiącego rozwiąza-
nie zadania regulacji predykcyjnej w wersji MPC. Zadanie to zostało przedstawione w ze-
stawieniu C. Do jego rozwiązania koniecznym jest ustalenie, oprócz modelu sterowanego
procesu, ciągów wag wskaz
nikowych (ywt,uwt ), a także horyzontu sterowania (M)
i predykcji (P).
Y
Y Przy pomocy m-funkcji:
 [y,u,ym]=cmpc(plant,model,ywt,uwt,M,P,tend,r,ulim,ylim,...)
możemy zasymulować pracę układu sterowania w którym, na bieżąco, cyklicznie rozwią-
zywane jest zadanie otrzymane w wyniku transformacji zadania sterowania predykcyjnego
(problem MPC z ograniczeniami stosowanych sterowań i uzyskiwanego wyjścia obiektu)
do postaci zadania programowania kwadratowego. Oprócz argumentów wymienionych już
w uprzednich rozważaniach przed wywołaniem programu danej funkcji należy ustalić czas
trwania procesu ( tend ), podać stałą lub zmieniającą się w czasie wartość zadaną ( r ),
a także (ewentualnie) dolne i górne wartości ( ulim ,  ulim ) stosowanych ograniczeń.
Pełna wersja tej m-funkcji zapewnia również możliwość symulacji zachowania się systemu
przy uwzględnieniu oddziaływania zadeklarowanych zakłóceń.
Wywołując m- funkcję:
 [y,u,ym]=mpcsim(plant,model,Kmpc,tend,r,usat,...)
możemy sprawdzić w wyniku symulacji przebiegów rezultat syntezy regulatora predyk-
cyjnego, otrzymanego za pośrednictwem  mpccon . Argumenty danej funkcji znane są
już objaśnień m- funkcji poprzednich  wprowadzana zmienna  usat przekazuje dolny
i górny poziomy nasycenia sygnału wyjściowego sterownika.
48
Należy podkreślić, że omawiany pakiet zawiera również (w blokach:  plotting and
matrix information i  utility functions ) szereg m-funkcji ułatwiających
wizualizację rezultatów badań symulacyjnych. Wymieńmy niektóre z nich:
plotall - zapewnia przedstawienie wszystkich otrzymanych przebiegów w układzie
sterowania predykcyjnego na jednym wykresie;
ploteach odwrotnie, zapewnia otrzymanie różnych przebiegów na oddzielnych wy-
kresach;
plotstep  przedstawia graficznie współczynniki modelu w  MPC step format ;
mpcstair  kreuje schodkową postać przebiegów sterowania, itp.
Wymienimy także grupy m-funkcji niewykorzystywane bezpośrednio w trakcie ćwiczeń:
m- funkcje bloku identyfikacji (  identification ) dotyczą głównie obiektów wie-
lowymiarowych;
m- funkcje  model building odgrywają rolę pomocniczą; wspomagają proces prze-
kształcania schematów blokowych;
m- funkcje ułatwiające analizę własności uzyskanego układu regulacji ( analysis ) -
w programie ćwiczenia z założenia stosowane są głównie charakterystyki czasowe
obiektów wyznaczane na podstawie modeli typu wejście-wyjście;
m- funkcje projektowania sterownika i symulacji układów uzyskiwanych przy zastosowa-
niu, sformułowanego w specyficznym zapisie (tzw. MPC mod format ), modelu
opartego na równaniach stanu ( controller design and simulation ).
49
UWAGI BIBLIOGRAFICZNE i LITERATURA.
Problemy identyfikacji obiektów dynamicznych w obecności zakłóceń stochastycznych
omawiane są w wielu podręcznikach. Znanym już z poprzednich zajęć jest [1]. Bardzo mia-
rodajną teoretycznie, a jednocześnie zawierającą wiele ważnych informacji pratycznych jest
monografia [2]. Podręczniki z zakresu sterowania adaptacyjnego i predykcyjnego napisane
w układzie bliskim reprezentowanemu w trakcie ćwiczeń pochodzą od autorów z Politechniki
Śląskiej [3], [4], [5]. Z książek w języku angielskim warto zapoznać się z przynajmniej jedną
z fundamentalnych prac K. J.Astroma, np. [6].
[1] P. de Larminat, Y. Thomas : Automatyka  układy liniowe. T.2, Identyfikacja. WNT,
Warszawa, 1983
[2]T. Soderstrom, P. Stoica: Identyfikacja systemów. Wyd. Naukowe PWN, Warszawa, 1997
[3] A. Niederliński, Systemy komputerowe automatyki przemysłowej, T.2, Zastosowania,
WNT, Warszawa, 1985
[4] A. Niederliński, J. Mościński , Z. Ogonowski: Regulacja adaptacyjna. PWN, Warszawa,
1995
[5] Jerzy Mościński, Antoni Niederliński, Pakiet Predal. Podręcznik użytkownika. Skrypt
uczelniany nr 1640 Politechniki Śląskiej. Gliwice,1991.
[6] K.J. Astrom, B. Wittenmark.: Adaptive Control, Addison-Wesley, 1989.
50