HW1: Init commit with assignment 1 source

.gitignore

@@ -0,0 +1,2 @@
+# Generic directories
+Inbox/

.gitmodules

@@ -0,0 +1,3 @@
+[submodule "Work 1/report/AUThReport"]
+	path = Work 1/report/AUThReport
+	url = ssh://git@git.hoo2.net:222/hoo2/AUThReport.git

Work 1/report/.gitignore

@@ -0,0 +1,6 @@
+# LaTeX auxiliary files
+*.aux
+*.log
+*.out
+*.synctex.gz
+

Work 1/report/AUThReport

@@ -0,0 +1 @@
+Subproject commit 74ec4b5f6c66382e5f1b6d2e6930897e4ed53ea6

Work 1/report/Work1_report.tex

@@ -0,0 +1,98 @@
+%
+% !TEX TS-program = xelatex
+% !TEX encoding   = UTF-8 Unicode
+% !TEX spellcheck = el-GR
+%
+% AUTH report template
+%
+% Requires compilation with pdfLaTeX or XeLaTeX
+%
+% authors:
+%   Χρήστος Χουτουρίδης ΑΕΜ 8997
+%   cchoutou@ece.auth.gr
+
+% Options:
+%
+% 1) mainlang=<language>
+%    Default:  english
+%    Set the default language of the document which affects hyphenations,
+%    localization (section, dates, etc...)
+%
+%    example: \documentclass[mainlang=greek]{AUThReport}
+%
+% 2) <language>
+%    Add hyphenation and typesetting support for other languages
+%    Currently supports: english, greek, german, frenc
+%
+%    example: \documentclass[english, greek]{AUThReport}
+%
+% 3) short: Requests a shorter title for the document
+%    Default: no short
+%
+%    example: \documentclass[short]{AUThReport}
+%
+\documentclass[a4paper, 11pt, mainlang=greek, english]{AUThReport/AUThReport}
+
+\CurrentDate{\today}
+
+% Greek report document setup suggestions
+%---------------------------------
+% \WorkGroup{Ομάδα Χ}
+
+\AuthorName{Χρήστος Χουτουρίδης}
+\AuthorAEM{8997}
+\AuthorMail{cchoutou@ece.auth.gr}
+
+%\CoAuthorName{Όνομα Επίθετο}
+%\CoAuthorAEM{1234}
+%\CoAuthorMail{xxx@ece.auth.gr}
+
+\DocTitle{Εργασία στα Ασαφή Συστήματα}
+\DocSubTitle{Έλεγχος γωνίας προσανατολισμού δορυφόρου με ασαφείς ελεγκτές}
+
+\Department{Τμήμα ΗΜΜΥ. Τομέας Ηλεκτρονικής}
+\ClassName{Ασαφή Συστήματα (Υπολογιστική Νοημοσύνη)}
+
+\InstructorName{Θεοχάρης Ιωάννης}
+\InstructorMail{theochar@ece.auth.gr}
+
+\CoInstructorName{Χαδουλός Χρήστος}
+\CoInstructorMail{christgc@auth.gr}
+
+
+\usepackage{float}
+\usepackage{minted}
+\usepackage{xcolor} %
+
+\setminted[cpp]{
+	fontsize=\small,
+	breaklines,
+	autogobble,
+	baselinestretch=1.1,
+	tabsize=2,
+	numbersep=8pt,
+	gobble=0
+}
+
+\newcommand{\repo}{https://git.hoo2.net/hoo2/PDS/src/branch/master/homework_3}
+
+\begin{document}
+
+% Request a title page or header
+\InsertTitle
+%\InsertTitle[img/background.png][0.8\textwidth][2cm]
+
+\section{Εισαγωγή}
+
+Η παρούσα εργασία αποτελεί
+
+
+\subsection{Παραδοτέα}
+Τα παραδοτέα της εργασίας αποτελούνται από:
+\begin{itemize}
+	\item Την παρούσα αναφορά.
+	\item Το \href{\hwThree}{σύνδεσμο με το αποθετήριο} που περιέχει τον κώδικα για την παραγωγή των εκτελέσιμων και της αναφοράς.
+\end{itemize}
+
+
+\end{document}

Work 1/source/mkFIS.m

@@ -0,0 +1,66 @@
+function fis = mkFIS()
+%mkFIS makes FIS
+
+
+
+% FIS control
+%  R2025 compatibility:
+%  - newfis->mamfis
+%  - defuzzMethod='cos' -> AggregationMethod="sum", DefuzzificationMethod="centroid" 
+fis = mamfis( ...
+    AndMethod="min", ...
+    OrMethod="max", ...
+    ImplicationMethod="min", ...
+    AggregationMethod="sum", ...
+    DefuzzificationMethod="centroid" );
+
+
+% Member functions
+mf.Names = ["NV"; "NL"; "NM"; "NS"; "ZR"; "PS"; "PM"; "PL"; "PV"];
+mf.Edges = [     ...
+    [-1.25 -1    -0.75]; ...  % NV
+    [-1    -0.75 -0.5 ]; ...  % NL 
+    [-0.75 -0.5  -0.25]; ...  % NM
+    [-0.5  -0.25  0   ]; ...  % NS
+    [-0.25  0     0.25]; ...  % ZR
+    [ 0     0.25  0.5 ]; ...  % PS
+    [ 0.25  0.5   0.75]; ...  % PM
+    [ 0.5   0.75  1   ]; ...  % PL
+    [ 0.75  1     1.25]; ...  % PV
+];
+
+% Setup input/output
+fis = addInput(fis, [-1 1], Name="e");
+fis = addInput(fis, [-1 1], Name="de");
+fis = addOutput(fis, [-1 1], Name="du");
+for k=1:length(mf.Names)
+    fis = addMF(fis, "e",  "trimf", mf.Edges(k, :), Name=mf.Names(k));
+    fis = addMF(fis, "de", "trimf", mf.Edges(k, :), Name=mf.Names(k));
+    fis = addMF(fis, "du", "trimf", mf.Edges(k, :), Name=mf.Names(k));
+end
+
+% Rules 9x9
+
+% rows = e (NV..PV), cols = de (NV..PV)
+ruleTable = [
+% de: NV   NL   NM   NS   ZR   PS   PM   PL   PV
+     "NV","NV","NV","NV","NV","NL","NM","NS","ZR";  % e = NV
+     "NV","NV","NV","NV","NL","NM","NS","ZR","PS";  % e = NL
+     "NV","NV","NV","NL","NM","NS","ZR","PS","PM";  % e = NM
+     "NV","NV","NL","NM","NS","ZR","PS","PM","PL";  % e = NS
+     "NV","NL","NM","NS","ZR","PS","PM","PL","PV";  % e = ZR
+     "NL","NM","NS","ZR","PS","PM","PL","PV","PV";  % e = PS
+     "NM","NS","ZR","PS","PM","PL","PV","PV","PV";  % e = PM
+     "NS","ZR","PS","PM","PL","PV","PV","PV","PV";  % e = PL
+     "ZR","PS","PM","PL","PV","PV","PV","PV","PV"   % e = PV
+];
+
+for i=1:length(mf.Names)
+    for j=1:length(mf.Names)
+        rule = "IF e IS "+mf.Names(i)+" AND de IS "+mf.Names(j)+" THEN du IS "+ruleTable(i,j);
+        fis = addRule(fis, rule);
+    end
+end
+
+
+end

Work 1/source/satFZ_PI.m

@@ -0,0 +1,45 @@
+function u = satFZ_PI(fis, e, Ke, Kde, Kdu, Ts)
+%SATFZ_PI Calculate the output of the controller
+% Must be called every Ts
+%
+% Output:
+%   u: Output
+
+% Inputs:
+%   fis: FIS to use
+%   e:   error
+%   Ke:  error gain
+%   Kde: error derivative gain = a*Ke
+%   Kdu: du gain
+%   Ts:  Sampling period
+
+persistent prev_e;
+persistent prev_u;
+
+if isempty(prev_e) % first call
+    prev_e = 0;
+    prev_u = 0;
+end
+
+% Saturation between min: m and max: M
+sat = @(x, m, M) max(m, min(M, x));
+
+
+% calculate error diff
+de  = (e - prev_e)/Ts;
+prev_e = e;
+
+% Normalize e and de
+en  = sat(Ke*e, -1, 1);
+den = sat(Kde*de, -1, 1);
+% evaluate FIS
+dun = evalfis(fis, [en den]);
+% De-normalize dun
+du  = Kdu * dun;
+
+
+% Update outputs
+u = du + prev_u;
+prev_u = u;
+
+end

Work 1/source/satelite_PI.m

@@ -0,0 +1,171 @@
+%% Satellite - Classical PI (Part A)
+%
+% Classical control script for Assignment 1 in Fuzzy systems
+% 
+% author:
+%   Christos Choutouridis ΑΕΜ 8997
+%   cchoutou@ece.auth.gr
+%
+% -------------------------------------------------------
+% Givens
+%                     10
+% Plant: Gp(s) = ------------
+%                 (s+1)(s+9)
+% 
+%                              s + c
+% PI: Gc(s) = Kp + Ki/s = Kp --------- , with c = Ki/Kp
+%                                s
+
+clear; clc; close all;
+
+s = tf('s');
+
+% Givens
+% -------------------------------------------------------
+Gp = 10/((s+1)*(s+9));          % given plant
+spec.Mp_max   = 0.10;           % 10% overshoot
+spec.Ts_max   = 2.0;            % 2sec settling time
+spec.SteadyStateValue = 1.0;    % step of 1
+
+% Configuration
+config.PI_zero = -1;
+config.stepresp_time = 5;
+
+% Control
+Gc = @(c)     (s + c)/s;            % No gain controller
+L0 = @(c)     0.1 * Gc(c) * Gp;     % No gain open loop
+Lk  = @(K, c) K * L0(c);            % open-loop with K [requested by the assignment]
+
+% Utilities to extract Kp, Ki from K [variable requested by the assignment]
+Kp = @(K)    0.1*K;
+Ki = @(K, c) c*Kp(K);
+
+% ---------- Brute force Root-locus exploration for fixed zero ----------
+c = -config.PI_zero;                 % zero location
+best = struct('cost',Inf);
+for K = linspace(0, 100, 400)
+    CL = feedback(Lk(K, c), 1, -1);  % closed loop unity feedback
+    info = stepinfo(CL);
+    if ~isfinite(info.SettlingTime) || isnan(info.Overshoot)
+        continue
+    end
+    % First check hard specs
+    if info.Overshoot/100 <= spec.Mp_max && ...
+       info.SettlingTime  <= spec.Ts_max
+        %
+        % Select the best CL based on a more appropriate cost function
+        % (ITAE-like): integral of t*|e(t)| approximated by weights
+        %
+        [y, t] = step(CL, 0:0.001:config.stepresp_time);
+        e = spec.SteadyStateValue - y;
+        cost = trapz(t, t.*abs(e));
+        if cost < best.cost
+            best.cost = cost;
+            best.K    = K;
+            best.info = info;
+            best.CL   = CL;
+            best.ess  = e(end);
+        end
+    end
+end
+
+fprintf('Best match: Kp=%.3g, Ki=%.3g\n', Kp(best.K), Ki(c, best.K));
+fprintf('Poles:\n');
+p = pole(best.CL);
+for k = 1:numel(p)
+    fprintf('  Pole %d: %.5f %+.5fi\n', k, real(p(k)), imag(p(k)));
+end
+fprintf('Step response:\n');
+fprintf('  RiseTime      : %.4g s\n', best.info.RiseTime);
+fprintf('  SettlingTime  : %.4g s\n', best.info.SettlingTime);
+fprintf('  Overshoot     : %.2f%%\n', best.info.Overshoot);
+fprintf('  Peak          : %.4g\n',   best.info.Peak);
+fprintf('  PeakTime      : %.4g s\n', best.info.PeakTime);
+fprintf('  SS error      : %.4g s\n', best.ess);
+
+f1 = figure( ...
+    'Color','w', ...
+    'Name',sprintf('Root Locus - c=%.3g', c), ...
+    'Position',[100 100 1200 800] ...
+);
+h = rlocusplot(L0(c));
+setoptions(h, 'XLim', [-12 2], 'YLim', [-8 8]);
+ax = findobj(h, 'Type','Axes');
+set(ax,'NextPlot','add');
+grid(ax,'on');
+
+% Overlays
+p_cl = pole(best.CL);
+hCL = plot(ax, real(p_cl), imag(p_cl), 'ro', 'MarkerSize', 8, ...
+           'MarkerFaceColor','r', 'DisplayName', ...
+           sprintf('CL poles: K_p=%.3g, K_i=%.3g', Kp(best.K), Ki(best.K, c)));
+
+infoRL = sprintf([ ...
+    'K_p = %.3g\n' ...
+    'K_i = %.3g\n' ...
+    'c  = %.3g\n' ...
+    'Poles:\n' ...
+    '  %.4f%+.4fi\n' ...
+    '  %.4f%+.4fi\n' ...
+    '  %.4f%+.4fi'], ...
+     Kp(best.K), Ki(c,best.K), c, ...
+     real(p_cl(1)), imag(p_cl(1)), ...
+     real(p_cl(2)), imag(p_cl(2)), ...
+     real(p_cl(3)), imag(p_cl(3)));
+
+p_ol = pole(L0(c));
+hOL = plot(ax, real(p_ol), imag(p_ol), 'kx', 'MarkerSize', 9, ...
+           'LineWidth', 1.4, 'DisplayName','OL poles');
+
+% z_ol = zero(L0(c));
+% hZ  = plot(ax, real(z_ol), imag(z_ol), 's', 'MarkerSize', 8, ...
+%            'MarkerFaceColor',[1 0.9 0.2], 'MarkerEdgeColor','k', ...
+%            'DisplayName','OL zero');
+
+% Legend - annotations
+lg = legend(ax, [hCL hOL], {'CL poles','OL poles'}, ...
+            'Location','best','AutoUpdate','off');
+
+annotation(f1, 'textbox', [0.75 0.05 0.4 0.3], ...
+    'String', infoRL, 'Interpreter','tex', ...
+    'FontName','monospaced', 'FontSize',12, ...
+    'FitBoxToText','on', 'BackgroundColor', 'w', 'EdgeColor', [1 1 1]);
+
+% titles
+% xlabel(ax,'Real Axis (s^{-1})'); ylabel(ax,'Imaginary Axis (s^{-1})');
+% title(ax, sprintf('Root Locus (c=%.3g)', c));
+drawnow;
+
+% save figure
+exportgraphics(f1, sprintf('Root Locus c%.3g.png', c), 'Resolution', 300);
+% exportgraphics(f1, sprintf('Root Locus c%.3g.pdf', c));
+
+
+f2 = figure( ...
+    'Color','w', ...
+    'Name',sprintf('Step responce - c=%.3g, Kp=%.3g, Ki=%.3g', c, Kp(best.K), Ki(best.K, c)), ...
+    'Position',[200 200 1200 800] ...
+);
+step(best.CL);
+xlabel('Time (s)');
+ylabel('Amplitude');
+title(sprintf('Step Response - c=%.3g, Kp=%.3g, Ki=%.3g', c, Kp(best.K), Ki(best.K, c)));
+grid on;
+
+infoStep = sprintf(['RiseTime   : %.4g s\n' ...
+                    'Settling   : %.4g s\n' ...
+                    'Overshoot  : %.2f %%\n' ...
+                    'Peak       : %.4g (at %.4g s)\n' ...
+                    'SS error   : %.3f'], ...
+                   best.info.RiseTime, best.info.SettlingTime, best.info.Overshoot, ...
+                   best.info.Peak, best.info.PeakTime, best.ess);
+
+annotation(f2, 'textbox', [0.7 0.05 0.4 0.25], ...
+    'String', infoStep, 'Interpreter','tex', ...
+    'FontName','monospaced', 'FontSize',12, ...
+    'FitBoxToText','on', 'BackgroundColor','w', 'EdgeColor',[1 1 1]);
+
+
+% save figure
+exportgraphics(f2, sprintf('Step Responce c%.3g_Kp%.3g_Ki%.3g.png', c, Kp(best.K), Ki(best.K, c)), 'Resolution', 300);
+% exportgraphics(f2, sprintf('Step Responce c%.3g_Kp%.3g_Ki%.3g.pdf', c, Kp(best.K), Ki(best.K, c)));

Work 1/source/scenario1_tune.m

@@ -0,0 +1,90 @@
+%% Satellite - Scenario 1 tune FZ-PI
+%
+% Assignment 1 in Fuzzy systems
+% 
+% author:
+%   Christos Choutouridis ΑΕΜ 8997
+%   cchoutou@ece.auth.gr
+%
+
+% 
+clear; clc; close all;
+
+s = tf('s');
+
+% Givens
+% -------------------------------------------------------
+Gp = 10/((s+1)*(s+9));          % given plant
+spec.Mp_max   = 0.07;           % 7% overshoot
+spec.Tr_max   = 0.6;            % 0.6sec rise time
+spec.SteadyStateValue = 1.0;    % step of 1
+
+Ts = 0.01;
+
+% Configuration
+tEnd = 5;
+Ke  = 1;
+fis = mkFIS();
+
+best = struct('cost',Inf);
+for Kde = 0.2:0.01:0.35
+    for Kdu = 0.3:0.01:0.5
+        % Convert Gp -> A,B,C,D (z)
+        Gz = c2d(Gp, Ts, 'zoh');
+        Gd = ss(Gz);
+        Ad=Gd.A; Bd=Gd.B; Cd=Gd.C; Dd=Gd.D;
+        
+        %figure('Name','FLC surface'); gensurf(fis); grid on;
+        
+        
+        t = 0:Ts:tEnd;
+        r = ones(size(t));
+        x = zeros(size(Ad,1),1);
+        y = 0;
+        
+        % logs
+        Y = zeros(size(t)); U = zeros(size(t)); E = zeros(size(t));
+        
+        clear satFZ_PI;  % Init satFZ_PI
+        for k = 1:numel(t)
+            % calculate error
+            e = r(k) - y;
+        
+            %evaluate controller
+            u = satFZ_PI(fis, e, Ke, Kde, Kdu, Ts);
+        
+            % plant update (ZOH)
+            x = Ad*x + Bd*u;
+            y = Cd*x + Dd*u;
+                
+            % keep logs
+            Y(k) = y; U(k) = u; E(k) = e;
+        end
+        info = stepinfo(Y, t, 1);
+        % First check hard specs
+        if info.Overshoot/100 <= spec.Mp_max && ...
+           info.RiseTime  <= spec.Tr_max
+            %
+            % Select the best CL based on a more appropriate cost function
+            % (ITAE-like): integral of t*|e(t)| approximated by weights
+            %
+            cost = trapz(t, t.*abs(E));
+            if cost < best.cost
+                best.cost = cost;
+                best.Kdu  = Kdu;
+                best.Kde  = Kde;
+                best.info = info;
+                best.ess  = E(end);
+            end
+        end
+    end
+end
+
+fprintf('Best match: Kde=%.3g, Kdu=%.3g\n', best.Kde, best.Kdu);
+fprintf('Step response:\n');
+fprintf('  RiseTime      : %.4g s\n', best.info.RiseTime);
+fprintf('  SettlingTime  : %.4g s\n', best.info.SettlingTime);
+fprintf('  Overshoot     : %.2f%%\n', best.info.Overshoot);
+fprintf('  Peak          : %.4g\n',   best.info.Peak);
+fprintf('  PeakTime      : %.4g s\n', best.info.PeakTime);
+fprintf('  SS error      : %.4g s\n', best.ess);

Work 1/source/scenario1a.m

@@ -0,0 +1,90 @@
+%% Satellite - Scenario 1a
+%
+% Assignment 1 in Fuzzy systems
+% 
+% author:
+%   Christos Choutouridis ΑΕΜ 8997
+%   cchoutou@ece.auth.gr
+%
+
+% 
+clear; clc; close all;
+
+s = tf('s');
+
+% Givens
+% -------------------------------------------------------
+Gp = 10/((s+1)*(s+9));          % given plant
+spec.Mp_max   = 0.07;           % 7% overshoot
+spec.Tr_max   = 0.6;            % 0.6sec rise time
+spec.SteadyStateValue = 1.0;    % step of 1
+
+Ts = 0.01;
+
+% Configuration
+tEnd = 5;
+Ke  = 1;
+Kde=0.29;   % Calculated with scenario1_tune.m
+Kdu=0.33;
+fis = mkFIS();
+
+Gz = c2d(Gp, Ts, 'zoh');
+Gd = ss(Gz);
+Ad=Gd.A; Bd=Gd.B; Cd=Gd.C; Dd=Gd.D;
+
+
+t = 0:Ts:tEnd;
+r = ones(size(t));
+x = zeros(size(Ad,1),1);
+y = 0;
+
+% logs
+Y = zeros(size(t)); U = zeros(size(t)); E = zeros(size(t));
+
+clear satFZ_PI;  % Init satFZ_PI
+for k = 1:numel(t)
+    % calculate error
+    e = r(k) - y;
+
+    %evaluate controller
+    u = satFZ_PI(fis, e, Ke, Kde, Kdu, Ts);
+
+    % plant update (ZOH)
+    x = Ad*x + Bd*u;
+    y = Cd*x + Dd*u;
+
+    % keep logs
+    Y(k) = y; U(k) = u; E(k) = e;
+end
+info = stepinfo(Y, t, 1);
+
+
+f1 = figure( ...
+    'Color','w', ...
+    'Name',sprintf('Step responce - Kde=%.3g, Kdu=%.3g', Kde, Kdu), ...
+    'Position',[200 200 1200 800] ...
+);
+plot(t, Y);
+xlabel('Time (s)');
+ylabel('Amplitude');
+title(sprintf('Step responce - Kde=%.3g, Kdu=%.3g', Kde, Kdu));
+grid on;
+
+infoStep = sprintf(['RiseTime   : %.4g s\n' ...
+                    'Settling   : %.4g s\n' ...
+                    'Overshoot  : %.2f %%\n' ...
+                    'Peak       : %.4g (at %.4g s)\n' ...
+                    'SS error   : %.3f'], ...
+                   info.RiseTime, info.SettlingTime, info.Overshoot, ...
+                   info.Peak, info.PeakTime, E(end));
+
+annotation(f1, 'textbox', [0.7 0.05 0.4 0.25], ...
+    'String', infoStep, 'Interpreter','tex', ...
+    'FontName','monospaced', 'FontSize',12, ...
+    'FitBoxToText','on', 'BackgroundColor','w', 'EdgeColor',[1 1 1]);
+
+
+% save figure
+exportgraphics(f1, sprintf('Step Responce_FZ.png'), 'Resolution', 300);
+% exportgraphics(f1, sprintf('Step Responce_FZ.pdf'));
+

Work 1/source/scenario1b_c.m

@@ -0,0 +1,25 @@
+%% Satellite - Scenario 1b, 1c
+%
+% Assignment 1 in Fuzzy systems
+% 
+% author:
+%   Christos Choutouridis ΑΕΜ 8997
+%   cchoutou@ece.auth.gr
+%
+
+% 
+clear; clc; close all;
+
+% load fuzzy controller
+flc = mkFIS();
+
+% test controller
+% e=NM  for ex: -0.41
+% de=ZR for ex: 0.03
+evalfis(flc, [0.41 0.03]);
+ruleview(flc);
+
+
+f1 = figure('Name','FLC surface'); gensurf(fis); grid on;
+exportgraphics(f1, sprintf('FLC surface.png'), 'Resolution', 300);
+% exportgraphics(f1, sprintf('FLC surface.pdf'));

Work 1/source/scenario2.m

@@ -0,0 +1,127 @@
+%% Satellite - Scenario 2
+%
+% Assignment 1 in Fuzzy systems
+% 
+% author:
+%   Christos Choutouridis ΑΕΜ 8997
+%   cchoutou@ece.auth.gr
+%
+
+% 
+clear; clc; close all;
+
+s = tf('s');
+
+% Givens
+% -------------------------------------------------------
+Gp = 10/((s+1)*(s+9));          % given plant
+spec.Mp_max   = 0.07;           % 7% overshoot
+spec.Tr_max   = 0.6;            % 0.6sec rise time
+spec.SteadyStateValue = 1.0;    % step of 1
+
+Ts = 0.01;
+
+% Configuration
+Ke  = 1;
+Kde=0.29;   % Calculated with scenario1_tune.m
+Kdu=0.33;
+fis = mkFIS();
+
+Gz = c2d(Gp, Ts, 'zoh');
+Gd = ss(Gz);
+Ad=Gd.A; Bd=Gd.B; Cd=Gd.C; Dd=Gd.D;
+
+
+% Control inputs
+tEnd = 16;
+t = 0:Ts:tEnd;
+
+% 2a
+r_step = zeros(size(t));
+r_step(t < 4)                 = 60;
+r_step(t >= 4 & t < 8)        = 20;
+r_step(t >= 8)                = 40;
+
+x = zeros(size(Ad,1),1);
+y = 0;
+% logs
+Y = zeros(size(t)); U = zeros(size(t)); E = zeros(size(t));
+
+clear satFZ_PI;  % Init satFZ_PI
+for k = 1:numel(t)
+    % calculate error
+    e = r_step(k)/60 - y;
+
+    %evaluate controller
+    u = satFZ_PI(fis, e, Ke, Kde, Kdu, Ts);
+
+    % plant update (ZOH)
+    x = Ad*x + Bd*u;
+    y = Cd*x + Dd*u;
+
+    % keep logs
+    Y(k) = y; U(k) = u; E(k) = e;
+end
+
+f1 = figure( ...
+    'Color','w', ...
+    'Name',sprintf('MultiStep responce'), ...
+    'Position',[200 200 1200 800] ...
+);
+plot(t, r_step, 'k--', 'LineWidth',1.2); hold on;
+plot(t, Y*60, 'b', 'LineWidth',1.2);
+grid on; xlabel('Time (s)'); ylabel('Amplitude');
+title('Scenario 2A: multi-step (60->20->40)');
+legend('r','y','Location','best');
+
+% save figure
+exportgraphics(f1, sprintf('MultiStep Responce_FZ.png'), 'Resolution', 300);
+% exportgraphics(f1, sprintf('Step Responce_FZ.pdf'));
+
+
+% 2b
+r_trap = zeros(size(t));
+ix = (t <= 5);
+r_trap(ix) = 60/5 * t(ix);
+ix = (t > 5 & t <= 8);
+r_trap(ix) = 60;
+ix = (t > 8 & t <= 16);
+r_trap(ix) = 60 * (1 - (t(ix) - 8)/(16 - 8));
+
+x = zeros(size(Ad,1),1);
+y = 0;
+% logs
+Y = zeros(size(t)); U = zeros(size(t)); E = zeros(size(t));
+
+clear satFZ_PI;  % Init satFZ_PI
+for k = 1:numel(t)
+    % calculate error
+    e = r_trap(k)/60 - y;
+
+    %evaluate controller
+    u = satFZ_PI(fis, e, Ke, Kde, Kdu, Ts);
+
+    % plant update (ZOH)
+    x = Ad*x + Bd*u;
+    y = Cd*x + Dd*u;
+
+    % keep logs
+    Y(k) = y; U(k) = u; E(k) = e;
+end
+info = stepinfo(Y, t, 1);
+
+
+f2 = figure( ...
+    'Color','w', ...
+    'Name',sprintf('Trapezoid responce'), ...
+    'Position',[200 200 1200 800] ...
+);
+plot(t, r_trap, 'k--', 'LineWidth',1.2); hold on;
+plot(t, Y*60, 'b', 'LineWidth',1.2);
+grid on; xlabel('Time (s)'); ylabel('Amplitude');
+title('Scenario 2B: trapezoid (0->60, hold, 60->0)');
+legend('r','y','Location','best');
+
+% save figure
+exportgraphics(f2, sprintf('Trapezoid Responce_FZ.png'), 'Resolution', 300);
+% exportgraphics(f2, sprintf('Step Responce_FZ.pdf'));