% =========================================================================
%   Praktikum MATLAB/Simulink I
%   Technische Universität Darmstadt
%   IAT, FG rtm (Prof. Dr.-Ing. U. Konigorski)
%
%   Copyright (c) 2017, FG Regelungstechnik und Mechatronik,
%                       Technische Universität Darmstadt
%   All rights reserved.
%
%   Redistribution and use in source and binary forms, with or withou
%   modification, are permitted provided that the following conditions are
%   met:
%
%   1. Redistributions of source code must retain the above copyright
%      notice, this list of conditions and the following disclaimer.
%
%   2. Redistributions in binary form must reproduce the above copyright
%      notice, this list of conditions and the following disclaimer in the
%      documentation and/or other materials provided with the distribution.
%
%
%   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
%   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
%   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
%   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
%   HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
%   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
%   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
%   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
%   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
%   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
%   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

% =========================================================================

% Versuch 4

clear all
clc


%% Versuchsdurchführung 4.1.1

% Übertragungsfunktion des Ventilators
sys = tf(0.006213, [8.989e-10 2.786e-5 4.71e-5]);

%% Versuchsdurchführung 4.1.2

% Bestimmung der Pole
p = pole(sys);
fprintf('Die Pole der Strecke sind %e und %e\n', p(1), p(2))

%% Versuchsdurchführung 4.1.3
% Nullstellen und Pole in der komplexen Ebene
x = (-3.5*10^4: 100: 0).';
y = zeros(size(x));
figure;
plot(x, y, ':', 'Color', 'black');
hold on;
plot(p(1), 0, 'x', p(2), 0, 'x', 'Color', 'black', 'linewidth', 2);
xlabel('Real Axis');
ylabel('Imaginary Axis');
xlim([-3.5e4 0]);
ylim([-1 1]);
title('Pole-Zero-Map');


%% Versuchsdurchführung 4.1.4
% stationäre Verstärkung
k = dcgain(sys);
fprintf('Die stationäre Verstärkung der Strecke ist %g\n', k);

%% Versuchsdurchführung 4.1.5
% Sprungantwort
[y1, x1] = step(sys);
x2 = (0:0.01:3.5).';
y2 = k*ones(size(x2));
figure();
plot(x2, y2,':', 'Color', 'black');
hold on;
plot(x1, y1, 'linewidth', 2)
xlabel('Time (sec)');
ylabel('Amplitude');
ylim([0 140]);
xlim([0 3.5]);
grid on;
title('Step Response');


%% Impulsantwort
[y1, x1] = impulse(sys);

figure();
plot(x1, y1, 'linewidth', 2)
xlabel('Time (sec)');
ylabel('Amplitude');
ylim([0 250]);
xlim([0 3.5]);
grid on;
title('Impulse Response');


%% Versuchsdurchführung 4.1.7
% Bodediagramm und Ortskurve

[mag, phase, wout] = bode(sys, 'w', {0.1, 1e6});
figure();
subplot(2, 1, 1);
semilogx(wout, 10*log(mag(:,:)), 'linewidth', 2);
grid on;
xlim([0.1 1e6]);
ylim([-150 50]);
ylabel('Magnitude (dB)');
set(gca(), 'XTickLabel', {''});
title('Bode Diagram');

subplot(2, 1, 2);
semilogx(wout, phase(:,:), 'linewidth', 2);
grid on;
xlim([0.1 1e6]);
ylim([-180 0]);
ylabel('Phase (deg)');
xlabel('Frequency (rad/sec)')


x1 = (-20:1:140).';
y1 = zeros(size(x1));
y2 = (-80:1:80).';
x2 = zeros(size(y2));
[x3, y3] = nyquist(sys);

figure();
plot(x1, y1, ':', 'Color', 'black');
hold on;
plot(x2, y2, ':', 'Color', 'black');
plot(-1, 0, 'x', 'Color', 'black')
plot(x3(:,:), y3(:,:), 'linewidth', 2);
plot(x3(:,:), -y3(:,:), 'linewidth', 2);
plot(65.4200, -65.95 , '<', 'Color', 'black', 'linewidth', 3);
plot(65.4200, 65.95, '>', 'Color', 'black', 'linewidth', 3);
xlim([-20 140]);
ylim([-80 80]);
xlabel('Real Axis');
ylabel('Imaginary Axis');
title('Nyquist Diagram');


[x3, y3] = nyquist(sys, 'w', {1e3, 1e9});

figure();
plot(x1, y1, ':', 'Color', 'black');
hold on;
plot(x2, y2, ':', 'Color', 'black');
plot(-1, 0, 'x', 'Color', 'black')
plot(x3(:,:), y3(:,:), 'linewidth', 2);
plot(x3(:,:), -y3(:,:), 'linewidth', 2);
plot(65.4200, -65.95 , '<', 'Color', 'black', 'linewidth', 3);
plot(65.4200, 65.95, '>', 'Color', 'black', 'linewidth', 3);
xlim([-0.04 0.003]);
ylim([-0.035 0.015]);
xlabel('Real Axis');
ylabel('Imaginary Axis');
title('Nyquist Diagram');


[Gm, phiR, Ard, WD] = margin(sys);
fprintf('omega_D=%g\nphi_R=%g\nA_rd=%g\n', WD, phiR, Ard)

%% Versuchsdurchführung 4.1.8

x1 = (0:0.1:16).';
y1 = zeros(size(x1));
u_sin=sin(x1);
u_cos=cos(x1);
[y_sin, x3]=lsim(sys,u_sin,x1);
[y_cos, x3]=lsim(sys,u_cos,x1);

figure;
plot(x1, y1, ':', 'Color', 'black');
hold on;
plot(x1, y_sin, 'linewidth', 2);
xlim([0 10]);
ylim([-150 150]);
xlabel('Time (sec)');
ylabel('Amplitude');
title('Linear Simulation Results');


%% Versuchsdurchführung 4.2.1.1

% Bodediagramm
figure('name', 'Bodediagramm')
bode(sys)
grid on

%% Versuchsdurchführung 4.2.1.3

% Reglerübertragungsfunktion
regler_P = 263.0268;

%% Versuchsdurchführung 4.2.1.4

% Übertragungsfunktion des geschlossenene Regelkreises
sys_ges = feedback(sys * regler_P, 1);

%Pole des geschlossenen Regelkreises
pg = pole(sys_ges);
fprintf('Pole des geschlossenene Regelkreises:\n');
pg

%% Versuchsdurchführung 4.2.1.5

% Reaktion des geschlossenen Regelkreises auf Einheitssprung
x1 = (0:0.5e-5:4e-4).';
y1 = ones(size(x1));
[y2, x2] = step(sys_ges, 4e-4);

figure();
plot(x2, y2, 'linewidth', 2);
legend('Geschlossener RK');
hold on;
plot(x1, y1, ':', 'Color', 'black');
xlabel('Time (sec)');
ylabel('Amplitude');
xlim([0 4e-4]);
ylim([0 1.4]);
title('Step Response');


%% Versuchsdurchführung 4.2.2.1

% PI-Reglerübertragungsfunktion
regler_PI_1 = tf([0.5917 1], [1 0]);

%% Versuchsdurchführung 4.2.2.2

% Bodediagramm zur Bestimmung der Durchtrittsfrequenz
figure('name', 'Bodediagramm zur Bestimmung der Durchtrittsfrequenz')
bode(sys * regler_PI_1)

%% Versuchsdurchführung 4.2.2.4

% Übertragungsfunktion des Reglers
regler_PI_1 = regler_PI_1 * 436.516

%% Versuchsdurchführung 4.2.2.5

% geschlossener Regelkreis
Gg = feedback(sys * regler_PI_1, 1);

% Reaktion des geschlossenen Regelkreises auf Einheitssprung
x1 = (0:0.5e-5:4e-4).';
y1 = ones(size(x1));
[y2, x2] = step(Gg, 4e-4);

figure();
plot(x2, y2, 'linewidth', 2);
legend('Geschlossener RK');
hold on;
plot(x1, y1, ':', 'Color', 'black');
xlabel('Time (sec)');
ylabel('Amplitude');
xlim([0 4e-4]);
ylim([0 1.4]);
title('Step Response');


%% Versuchsdurchführung 4.2.2.6

% Pole des geschlossenen Regelkreises
fprintf('Pole des geschlossenen PI-Regelkreises:')
pg = pole(Gg)

%% Versuchsdurchführunge 4.2.2.7

% Kompensieren des zweiten Pols
regler_PI_12 = tf([1/(3.1e4) 1],[1 0]);
Gg = feedback(sys * regler_PI_12 * 0.025, 1);


x1 = (0:0.1:7).';
y1 = ones(size(x1));
[y2, x2] = step(Gg, 7);

figure('name','PI-Regler: Sprungantwort bei Kompensierung des zweiten Pols')
plot(x2, y2, 'linewidth', 2);
legend('Geschlossener RK');
hold on;
plot(x1, y1, ':', 'Color', 'black');
xlabel('Time (sec)');
ylabel('Amplitude');
xlim([0 7]);
ylim([0 1.4]);
title('Step Response');



%% Versuchsdurchführung 4.3.3

% PI-Reglerübertragungsfunktion für PI Regler bei Reglersynthese nach Betragsoptimum
r0 = 234.9;   % Parameter r0
r1 = 138.9;   % Parameter r1
regler_PI_2 = tf([r1 r0], [2 0]);

%% Versuchsdurchführung 4.3.4

% Reaktion des geschlossenen Regelkreises bei Einheitssprung
sys_ges = feedback(sys * regler_PI_2, 1);

x1 = (0:0.5e-5:4e-4).';
y1 = ones(size(x1));
[y2,x2] = step(sys_ges, 4e-4);

figure('name','Betragsoptimum: Sprungantwort');
plot(x2, y2, 'linewidth', 2);
hold on;
plot(x1, y1, ':', 'Color', 'black');
xlabel('Time (sec)');
ylabel('Amplitude');
xlim([0 4e-4]);
ylim([0 1.4]);
title('Step Response');



% Verdopplung von r0
regler_PI_22 = tf([r1 2*r0], [2 0]);
sys_ges_mod = feedback(sys * regler_PI_22, 1);

x1 = (0:0.5e-5:4e-4).';
y1 = ones(size(x1));
[y2, x2] = step(sys_ges_mod, 4e-4);

figure('name', 'Betragsoptimum: Sprungantwort bei Verdopplung von r0')
plot(x2, y2, 'linewidth', 2);
hold on;
plot(x1, y1, ':', 'Color', 'black');
xlabel('Time (sec)');
ylabel('Amplitude');
xlim([0 4e-4]);
ylim([0 1.4]);
title('Step Response');

disp('Pole des alten Systems:')
pole(sys_ges)
disp('Pole des neuen Systems:')
pole(sys_ges_mod)

%% Versuchsdurchführung 4.4.3

% reales Messglied
mess = tf(1, [0.00001 1]);


figure();
sys_ges2 = feedback(sys * regler_P, mess);
step(sys_ges2);
legend('P-Regler')

figure();
sys_ges2 = feedback(sys * regler_PI_1, mess);
step(sys_ges2);
legend('PI-Regler')

figure();
sys_ges2 = feedback(sys * regler_PI_2, mess);
step(sys_ges2);
legend('PI-Regler Betragsopt')

