%% Návrh diskrétního regulátoru - cvičení
%% Zadání
% 1) Pro zadanou soustavu navrhněte diskrétní korekční člen (FR1) splňující podmínku nulové ustálené odchylky a konečného regulačního děje v okamžicích vzorkování i mimo ně, pro jednotkový skok řízení.
% 2) Navrhněte druhý korekční člen (FR2) tak, aby první výstupní hodnota Fw byla Fw(1)=0.5.
% 3) Navrhněte korekční člen (FR3) splňující podmínku nulové ustálené odchylky a konečného regulačního děje v okamžicích vzorkování i mimo ně, pro jednotkový skok poruchy.
% 4) Srovnejte dosažené řešení s výsledky dosaženými řešením pomocí PID a PSD regulátoru.
% Předpokládejte standardní regulační zpětnovazební zapojení se vstupem poruchy na vstupu soustavy.
% Návrh proveďte na základě dodaného řešení pro jednoduchou soustavu. Prostudujte si řešení pro jednoduchou soustavu a převeďte na řešení soustavy zadané.
% Jednoduchá soustava pro ukázku řešení:    F\left(s\right)=\frac{1}{s+1}
% Zadaná soustava :                                        F\left(s\right)=\frac{10s+1}{\left(30s+1\right)\left(8s+1\right)\left(3s+1\right)}
%
%% Řešení
% Volba periody vzorkování: čím bude menší, tím více se bude výsledek podobat spojitému průběhu. Malá perioda vzorkování ovšem vyžaduje dražší HW (rychlejší D/A a A/D převodníky, rychlejší procesor -> dražší řešení). Menší perioda vzorkování vede též k většímu akčnímu zásahu (viz. níže). Perioda vzorkování by měla být alespoň taková, aby se soustava chovala jako dolnofrekvenční propust vůči vstupnímu signálu, aby vyfiltrovala vstupní (schodovitý) signál (to uskutečňuje i největší časová konstanta. Čitatele nebo jmenovatele?).
%% Úvodní nastavení

close all
clear all
clc

FS = tf(1,[1 1])    % zadaná soustava
TV = 0.1;           % perioda vzorkování
TMAXGRF0 = 8;       % maximální čas pro zobrazování grafů otevřené smyčky
TSM = 3;            % doba simulace - přenosy uzavřené smyčky
model_name = 'cv_diskr_reg_21a.slx'; % jméno příslušného modelu SIMILINK

% V dodaném souboru zapojte modely obvodů pro testování regulátorů navržených pro 
% řízení a poruchu.

% převod na diskrétní soustavu (stejné hodnoty v okamžicích vzorkování 
% jako spojitá => vhodná pro návrh diskrétního regulátoru).
% Tvarovač nultého řádu zajistí po částech spojitý signál na vstupu
% soustavy (pro většinu soustav lepší než dirakovy pulzy, energie se
% rozloží do celé periody => lepší i pro akční člen).
FC = c2d (FS,TV,'zoh')

% srovnání odezev (pro soustavy s integrátorem je lepší vykreslení pomocí impulse)
figure(1)
step(FS,TMAXGRF0)
title('Srovnani odezvy spojiteho prubehu a diskretniho ekvivalentu')
hold on
step(FC)

% Při převodu do diskrétního tvaru se v přenosu objeví nula, jejíž pozice se mění s periodou vzorkování, a může se objevit i mimo jednotkovou kružnici. To může znepříjemnit návrh. Většinou při zmenšování periody vzorkování dojde ke vzdálení nuly od jednotkové kružnice.
% U impulsové charakteristiky je problém s definicí (či spíše filozofií použití) dirakova (im)pulsu. Když v následujícím srovnáte impulsové charakteristiky, budou tvarově stejné. Ale pokud sečteme všechny členy impulsové charakteristiky, nedojdeme k hodnotě jedna (jak by se očekávalo - přechodová charakteristika končí v jedné a to je suma členů impulsové charakteristiky), ale výsledkem bude hodnota 1/TV - abychom dosáhli hodnoty jedna, musíme průběh vynásobit periodou vzorkování. Důvodem je to, že MATLAB negeneruje impulsy, ale charakteristiku, která má stejnou plochu (energii) jako ekvivalentní spojitá (proto jsou tvary stejné, ale pulsy v součtu nedávají hodnotu jedna). V případě, že potřebujeme získat "naši/skutečnou/matematickou" impulsní charakteristiku, musíme hodnoty získané z příkazu impulse vynásobit periodou vzorkování ("energie" nebude v celé periodě, ale pouze v "impulsu").


figure(11)
impulse(FS,TMAXGRF0)
title('Srovnani spojiteho a diskretniho impulsu ')
hold on
impulse(FC)
soucet=sum( impulse(FC))

% zjištění čitatelového polynomu přenosu FC
[CIT_FC, JMEN_FC] = tfdata(FC,'v'); % rozklad na citatel a jmenovatel
SUM_CIT = sum(CIT_FC); % získání hodnoty pro normalizaci koeficientů na součet jedna ("příjemnější" výpočty pro lidi)
CIT_FC_NORM = CIT_FC / SUM_CIT;
CIT_FC
CIT_FC_NORM
% vezměte jeden z výše uvedených čitatelových polynomů a z něj vypočítěte 
% (na papíře) přenos FW (výsledek by měl být stejný pro oba polynomy)

%% Podmínky na tvorbu přenosu řízení Fw
% Čitatel diskretizovaného přenosu se nazývá P (lepší získat ve tvaru se zápornými mocninami, MATLAB je v kladných) :   \textrm{Fc}\left(z^{-1} \right)=\frac{P\left(z^{-1} \right)}{\ldotp \ldotp \ldotp }
% Možné podmínky pro výpočet regulátoru:
% a) Stabilita Fw z hodnot čitatele FC: přenos Fw musí obsahovat nuly FC ležící mimo jednotkovou kružnici.
% b) Stabilita Fw z hodnot jmenovatele FC: přenos 1-Fw musí obsahovat póly FC ležící mimo jednotkovou kružnici.
% Pozn.: stabilní spojitá soustava má stabilní diskrétní ekvivalent (nestabilní spojitá má nestabilní diskrétní). I když spojitá soustava nemá nuly, diskrétní ekvivalent nuly mít může (poloha nuly je funkcí periody vzorkování, čím je Tv menší, tím je nula dále od jednotkové kružnice a tím obtížnější může být návrh).
% c) Podmínka pro konečný regulační děj mimo okamžiky vzorkování: přenos Fw obsahuje čitatelový polynom P přenosu Fc. Tato podmínka v sobě obsahuje i podmínku stability (bod a).
% d) Podmínka nulové ustálené odchylky pro vstupní signál: přenos 1-Fw musí obsahovat člen (1-z^-1)^i, kde i je řád vstupního signálu (i=1 pro jednotkový skok ...)
% Jelikož výše uvedené podmínky musí platit současně, je nutné pro řešení Fw podmínky zapsat společně a přidat polynomy, které zajistí současné splnění podmínek (umožní vyřešit rovnici) - polynomy R a M. Výsledné Fw bude mít tolik členů, kolik je nejjednodušší řešení rovnice. Dalšími členy lze upravit průběh Fw, ale bude mít více členů (bude delší). Z podmínky realizovatelnosti a tvaru zapojení je patrné, že pokud regulátor "nevěští", budou mít přenosy Fw a Fu stejné zpoždění (bude jim chybět tolik prvních členů) jako FC (u regulátoru předpokládáme, že bude mít odezvu již v nultém kroku).
% Sestavení rovnice ve tvaru:
% základní rovnice = hledané koeficienty Fw = podmínky pro čitatele = podmínky pro jmenovatele.
% Rozpis na koeficienty:
% 1-F_w =1-\left(f_0 +f_1 z^{-1} +f_2 z^{-2} +\ldotp \ldotp \ldotp \right)=1-\left(P\left(z^{-1} \right)R\left(z^{-1} \right)\right)=\left(\textrm{kořeny}\;\textrm{mimo}\;\textrm{jednotkovou}\;\textrm{kružnici}\right){\left(1-z^{-1} \right)}^i M\left(z^{-1} \right)
% V předchozí rovnici hledáme vhodné koeficienty polynomů R a M, tak aby vyhovovala poslední rovnost. Následně pomocí druhé rovnosti najdeme hodnoty pro koeficienty f_i. Koeficienty f_i jsou koeficienty přenosu řízení F_w. V případě, že máme požadavek na hodnotu přechodové charakteristiky (například hodnotu v prvním kroku řízení t=TV), tak musíme pro řešení využít poslední i předposlední rovnost (výše uvedené rovnice) společně, tak, že položíme f_1 =\textrm{"požadovaná}\;\textrm{hodnota}\;v\;\textrm{daném}\;\textrm{kroku"}.
% Pro naše řešení můžeme rozepsat na jednotlivé koeficienty:
% 1-F_w =1-\left(f_0 +f_1 z^{-1} +f_2 z^{-2} +\ldotp \ldotp \ldotp \right)=1-\left(\;1z^{-1} \;\ldotp \;\left({r_0 \;+r}_1 \;z^{-1} \;{+r}_2 \;z^{-2} +\ldotp \ldotp \ldotp \right)\right)={\left(1-z^{-1} \right)}^1 \left({m_0 \;+m}_1 \;z^{-1} \;{+m}_2 \;z^{-2} +\ldotp \ldotp \ldotp \right)
% Polynomy jsou si rovny, pokud se rovnají koeficienty u členů se stejnými mocninami:
% z^0:   1 - f_0 = 1  = m_0
% z^{-1}: -f_1= - r_0 = -m_0 \;     protože tato soustava rovnic stačí k vyřešení všech členů, dále již počítat nemusíme (tj nepřidáváme více členů polynomu M a R než je nutné).
% Řešením je:
% {F_w =f_1 z^{-1} \;=\;1\ldotp z}^{-1}
% Po vyřešení přenosu řízení se vypočítá regulátor po úpravě vzorce pro přenos řízení (vzorkování mezi všemi přenosy)
% F_w =\frac{\textrm{FC}\ldotp \textrm{FR}}{1+\textrm{FC}\ldotp \textrm{FR}}      =>      FR = ...
% Pozn.: v případě, že nebudeme (mít možnost) přenášet FR binárně, měli bychom se při zaokrouhlování jeho koeficientů snažit postupovat tak, abychom pokud možno dodrželi výše uvedené podmínky.
%


%%

% SIMULINK používá přenosy: FS a FRM pro simulaci; 
% pro stanovení konečného času simulace používá proměnnou TSM
% Aby SIMULINKový model znal proměnné z tohoto skriptu, je nutné skript spustit
% a potom použít simulinkovský model (teprve potom nahradí otazníky v blocích
% skutečné názvy).



disp('---------- kratky prubeh navrh rizeni -----------------')

FW1 = tf(1,[1 0],TV)

FR1 = FW1 / (FC *(1-FW1))
zpk(FR1)
% z přenosu FR je vidět, že se nezkrátily hodnoty, ale také výskyt členu
% (1-z^-1) ve jmenovateli, který zajišťuje nulovou ustálenou odchylku

% kontrola - zda po výpočtech "dokola" vyjde stejné FW jako na začátku
FWK1 = FR1*FC / (1+FR1*FC)
% (asi to vyjde jinak=>) můžeme se na to podívat z pohledu kořenů
zpk(FWK1)
% pokud je vidět, že se jedná o chyby na vyšších řádech -> zjednodušit
% přenos pomocí minreal - mělo by již vyjít
FWM1=minreal(FWK1)

% nyní můžeme vypočítat a zobrazit důležité přenosy uzavřené smyčky -
% řízení (máme);  poruchy; akčního zásahu
FUR1 = FC  / (1 + FC * FR1);
FXR1 = FR1 / (1 + FC * FR1);

figure(2);
step(FW1,TSM);
hold on
step(FUR1);
title({'Odezva rizeni a poruchy na jednotkovy skok' ; 'spojity a diskretizovany prubeh'})
% je vidět, že i když je řízení konečné v jednom kroku, porucha odeznívá
% delší dobu

% spojitý průběh řízení ze simulinku
FRM=FR1;
out = sim(model_name);
plot(out.ScopePorucha.time,out.ScopePorucha.signals.values  )
plot(out.ScopeRizeni.time  ,out.ScopeRizeni.signals(2).values )

figure(3)
step(FXR1,TSM)
title('Odezva akcniho zasahu pro rizeni na jednotkovy skok')
% zde je vidět jak je maximální hodnota akčního zásahu řádově větší než je
% hodnota "běžná" pro udržení ustálené odchylky

% Pozn.: Soustava je stále spojitá, i když pro výpočty používáme její diskretizovanou podobu, protože v uzavřené smyčce se vyskytují pouze hodnoty v okamžicích vzorkování. V grafech jsou zobrazeny oba průběhy, diskretizovaný i spojitý, který získáme simulací v SIMULINKu.
% Řešení s doplňkovými podmínkami:
% Řešení můžeme rozšířit dalšími požadavky, například, aby nebyl skok řízení tak prudký, ale byl pozvolnější, čehož dosáhneme rozložením přechodové charakteristiky do více kroků. Např.
% {F_w =}^{\;} \;{f_1 \;z}^{-1} {+f_2 \;z}^{-2} ={0,5\;z}^{-1} {+\;0,5\;z}^{-2} =\left(z^{-1} {+z}^{-2} \right)/\;2
% Z tohoto přenosu řízení se opět vypočítá regulátor a následně přenos řízení, poruchy a akční veličiny.
% Ke stejnému výsledku bychom došli i v případě dodatečné podmínky f_1 =0,5
% 1-F_w =1-\left(f_0 +f_1 \;z^{-1} +f_2 z^{-2} +\ldotp \ldotp \ldotp \right)=1-\left(f_0 +0,5\;z^{-1} +f_2 z^{-2} +\ldotp \ldotp \ldotp \right)=1-\left(\;1z^{-1} \;\ldotp \;\left({r_0 \;+r}_1 \;z^{-1} \;{+r}_2 \;z^{-2} +\ldotp \ldotp \ldotp \right)\right)={\left(1-z^{-1} \right)}^1 \left({m_0 \;+m}_1 \;z^{-1} \;{+m}_2 \;z^{-2} +\ldotp \ldotp \ldotp \right)
% z^0:   1 - f_0 = 1  = m_0
% z^{-1}: -f_1=-0,5= - r_0 = -m_0 \;+m_1 \;
% z^{-2}: -f_2= - r_1 = -m_1 \;     protože tato soustava rovnic (tento počet koeficientů m a r) stačí k vyřešení všech členů, dále již nemusíme rozšiřovat rovnice o další členy.
% F_w =\frac{\textrm{FC}\ldotp \textrm{FR}}{1+\textrm{FC}\ldotp \textrm{FR}}       FR = ...
% Výsledné průběhy zobrazíme do původních grafů
%


%%
disp('----------- delsi prubeh navrh rizeni ----------------')
FW2 = tf([1 1]/2,[1 0 0],TV)

FR2 = FW2 / (FC *(1-FW2))
zpk(FR2)
% z přenosu FR je vidět opět výskyt členu
% (1-z^-1) ve jmenovateli, který zajišťuje nulovou ustálenou odchylku
% přenos je složitější oproti původnímu s kratší Fw

% kontrola - zda po výpočtech "dokola" vyjde stejné FW jako na začátku
FWK2 = FR2*FC / (1+FR2*FC)
% (asi to vyjde jinak=>) můžeme se na to podívat z pohledu kořenů
zpk(FWK2)
% pokud je vidět, že se jedná o chyby na vyšších řádech -> zjednodušit
% pomocí minreal - mělo by již vyjít
FWMa=minreal(FWK2,0.0000001)

% nyní můžeme vypočítat a zobrazit důležité přenosy uzavřené smyčky -
% řízení (máme);  poruchy; akčního zásahu
FUR2 = FC /(1+FC*FR2);
FXR2 = FR2/(1+FC*FR2);

figure(2);
step(FW2,TSM);
title('Odezva rizeni a poruchy na jednotkovy skok - oba prenosy FR')
hold on
step(FUR2);
% je vidět, že i když je řízení konečné ve dvou krocích, porucha odeznívá
% delší dobu

% spojitý průběh získaný ze simulinku
FRM = FR2;
out = sim(model_name);
plot(out.ScopePorucha.time,out.ScopePorucha.signals.values  )
plot(out.ScopeRizeni.time  ,out.ScopeRizeni.signals(2).values )

figure(3)
hold on
step(FXR2,TSM)
title('Odezva akcniho zasahu na rizeni na jednotkovy skok - oba prenosy FR')
% zde je vidět jak je maximální hodnota akčního zásahu se prodloužením
% přechodného děje snížila



%% Řešení regulátoru na skok poruchy
% Jak jsme viděli v předchozím (protože nebyl integrační člen v regulátoru), nedošlo ke kvalitní regulaci poruchy. Můžeme tedy provádět návrh i s podmínkami na regulaci poruchy. V případě současných požadavků na řízení i poruchu by bylo nutné navrhnout regulátory dva (v přímé i zpětné větvi).
% Pro návrh na poruchu platí podobná pravidla jako pro řízení, pouze tam, kde byl člen 1-Fw je pouze člen Fu. To značně zjednoduší řešení (pokud nemáme další podmínky/požadavky).
% F_u =\left(f_0 +f_1 z^{-1} +f_2 z^{-2} +\ldotp \ldotp \ldotp \right)=\left(P\left(z^{-1} \right)R\left(z^{-1} \right)\right)=\left(\textrm{kořeny}\;\textrm{mimo}\;\textrm{jednotkovou}\;\textrm{kružnici}\right){\left(1-z^{-1} \right)}^i M\left(z^{-1} \right)
% Pozn.: Čitatelový polynom P zde používáme nenormalizovaný. Zatímco u návrhu na řízení je výhodou normalizace polynomu P, u návrhu na poruchu je nutné použít polynom P nenormalizovaný (díky podmínce, že první člen přechodové charakteristiky musí být stejný jako první člen P).
% Z uvedeného je vidět, že pro stabilní soustavu je řešením to, že zvolíme neznámý polynom tak, aby se rovnal známému na druhé straně rovnice
% tj.: M = P                     a                    R\left(z^{-1} \right)={\left(1-z^{-1} \right)}^i
% Pro uvedený příklad je řešením {F_u =1-z}^{-1}
% Po vyřešení přenosu řízení se vypočítá regulátor po úpravě vzorce pro přenos řízení (vzorkování mezi všemi přenosy)
% F_u =\frac{\textrm{FC}}{1+\textrm{FC}\ldotp \textrm{FR}}       FR = ...
%


%%
disp('---------- prubeh navrh porucha -----------------')

%[P, Q] = tfdata(FC,'v');
P = CIT_FC;

% při výpočtu se první člen FU odečítá s prvním členem polynomu P(Z).
% Pokud tyto členy nejsou stejné (např. jeden je zaokrouhlen), může se stát,
% že nedojde k jejich správnému odečtení (tj. pro odečtení zůstane část,
% kterou jsme zaokrouhlili) a regulátor "věští".
% Proto je pro výpočet v MATLABu nutné získat přesnou hodnotu (nebo
% vypočítat ručně; ...)
FU = tf([1 -1],[1 0],TV)*tf(P,[1 0],TV);

FR3 = 1/FU -  1 / FC
zpk(FR3)
% z přenosu FR je vidět výskyt členu (1-z^-1) ve jmenovateli, který zajišťuje nulovou ustálenou odchylku

% kontrola - zda po výpočtech "dokola" vyjde stejné FU jako na začátku
FUK = FC / (1+FR3*FC)
% (asi to vyjde jinak=>) můžeme se na to podívat z pohledu kořenů
zpk(FUK)
% % pokud je vidět, že se jedná o chyby na vyšších řádech -> zjednodušit
% % pomocí minreal - mělo by již vyjít
FUM=minreal(FUK)
%
% nyní můžeme vypočítat a zobrazit důležité přenosy uzavřené smyčky -
% řízení ;  poruchy(máme); akčního zásahu
FW3 = FR3 * FC/(1+FC*FR3);
FXR3= FR3     /(1+FC*FR3);
%
figure(4);
step(FW3,TSM);
title('Odezva rizeni a poruchy na jednotkovy skok')
hold on
step(FU);
% je vidět, že porucha je konečná v jednom kroku, naopak řízení odeznívá
% déle a má i horší průběh (díky překmitu), ale díky sumační složce v regulátoru
% a jednoduchému přenosu je i odezva řízení kvalitní

FRM = FR3;
out = sim(model_name);
plot(out.ScopePorucha.time,out.ScopePorucha.signals.values  )
plot(out.ScopeRizeni.time  ,out.ScopeRizeni.signals(2).values )

figure(5);
step(FXR3,TSM);
title('Odezva akcniho zasahu na jednotkovy skok rizeni')
% zde je vidět , že i zde je maximální hodnota akčního zásahu řádově větší než je
% hodnota "běžná" pro udržení ustálené odchylky

% Pokud ma regulator nulu mimo jednotkovou kruznici, ktera vytvari bipol,
% stanovte a použijte vhodnější periodu vzorkovani


% Často se stane, že má regulátor pól mimo jednotkovou kružnici. Pokud se jedná o bipól (c nulou v soustavě), potom je nutné se této nuly zbavit. To lze často provést změnou periody vzorkování. U některých soustav je to ovšem možné za cenu velice velkých period vzorkování - je ovšem lepší mít soustavu pomalou a stabilní, než nestabilní.