% Page108：已知闭环传递函数，求单位阶跃响应和单位斜坡响应
clear;
clc;

num = [1];
den = [1, 0.4, 1];
% 响应时间
t = [0:0.1:10];
% u 为单位斜坡输入
u = t;
% 单位阶跃响应
y = step(num, den, t);
% 单位斜坡响应
y1 = lsim(num, den, u, t);
% 将两条响应曲线绘制在同一个图上
plot(t, y, 'b-', t, y1, 'r:')
% 添加栅格
grid
% 添加横纵坐标轴
xlabel('Time [ses] t'); ylabel('y');
% 添加图标题
title('单位阶跃和单位斜坡输入响应曲线');
% 添加文字标注
legend('单位阶跃响应曲线', '单位斜坡响应曲线');

% Page109：已知开环传函，系统输入信号为三角波，求响应
clear;
clc;

numg = [1, 2];
deng = [1, 10, 1];
% 建立单位负反馈传递函数
[num, den] = cloop(numg, deng, -1);
% 产生三角波
v1 = [0:0.1:1];
v2 = [0.9:-0.1:-1];
v3 = [-0.9:0.1:0];
u = [v1, v2, v3];
% 仿真时间
t = [0:0.1:4];
% 求取系统在三角波输入下的响应
[y, x] = lsim(num, den, u, t);
% 将输入的三角波和输出响应曲线绘制在同一个图上
plot(t, y, t, u);
grid
% 添加横纵坐标轴
xlabel('Time [ses] t'); ylabel('y');
% 添加图标题
title('三角波输入和输入下的响应曲线');
% 添加文字标注
legend('响应曲线', '三角波输入');

% Page112：使用游动鼠标法和编程法求取系统的性能指标
clear;
clc;

% 建立零极点模型
G = zpk([ ], [-1+3i, -1-3i], 3);
% 求取阶跃响应
step(G);

% 以上为游动鼠标法
% 以下为编程法

% 建立零极点模型
G = zpk([ ], [-1+3i, -1-3i], 3);
% 计算最大峰值时间和它对应的超调量
C = dcgain(G);
% 求取阶跃响应
[y, t] = step(G);
plot(t,y)
grid
% 取得最大峰值时间
[Y, k] = max(y);
timetopeak = t(k);
% 计算超调量
percentovershoot = 100 * (Y-C)/C;
% 计算上升时间
n = 1;
while y(n) < C
    n = n+1;
end
risetime = t(n);
% 计算稳态响应时间
i = length(t);
while (y(i)>0.98*C)&(y(i)<1.02*C)
    i = i-1;
end
settlingtime = t(i);

% Page114：已知开环传函，计算阻尼比不同时单位负反馈系统的单位阶跃响应曲线
clear;
clc;

% 固有频率
wn = 1;
% 7 个不同的阻尼比取值
sigma = [0, 0.2, 0.4, 0.6, 0.9, 1.2, 1.5];
num = wn*wn;
% 将 t 在 0 到 20 之间均等分成 200 份
t = linspace(0, 20, 200);
for j = 1:7
    % 求开环传函的分母
    den = conv([1,0], [1,2*wn*sigma(j)]);
    % 建立开环传递函数
    s1 = tf(num, den);
    % 建立单位负反馈系统的传递函数
    sys = feedback(s1, 1);
    % 求取单位阶跃响应
    y(:, j) = step(sys, t);
end
% 在同一图上绘制单位阶跃响应曲线并添加栅格
plot(t, y(:, 1:7));grid
% 添加图标题
title('典型二阶系统取不同阻尼比时的单位阶跃响应');
% 放置 sigma 取不同值的文字注释
gtext('sigma=0');
gtext('sigma=0.2');
gtext('sigma=0.4');
gtext('sigma=0.6');
gtext('sigma=0.9');
gtext('sigma=1.2');
gtext('sigma=1.5');

% Page116：已知开环传函，求k不同时单位负反馈系统的单位阶跃响应曲线
clear;
clc;

% 时间常数
T = 1;
% 6 个不同的开环增益取值
k = [0.1, 0.2, 0.5, 0.8, 1.0, 2.4];
% 将时间 t 在 0 到 20 之间分成 200 份
t = linspace(0, 20, 200);
% 开环传递函数的分子、分母表达式
num = 1;
den = conv([1, 0], [T 1]);
for j = 1:6
    % 建立开环传递函数
    s1 = tf(num*k(j), den);
    % 建立闭环传递函数
    sys = feedback(s1, 1);
    % 求取单位阶跃响应
    y(:, j) = step(sys, t);
end
% 在同一图上绘制单位阶跃响应曲线并添加栅格
plot(t, y(:, 1:6));
grid
% 添加图标题
title('典型二阶系统取不同开环增益时的单位阶跃响应')
% 放置 k 取不同值的文字注释
gtext('k=0.1');
gtext('k=0.2');
gtext('k=0.5');
gtext('k=0.8');
gtext('k=1.0');
gtext('k=2.4');

% Page117：已知二阶系统，绘制tau不同值时，单位负反馈系统的单位阶跃曲线
clear;
clc;

% 时间常数
T = 1;
% 开环增益 K 的值
K = 1;
% 6 个不同的微分反馈系数
tau = [0, 0.05, 0.2, 0.5, 1.0, 2.4];
% 将 t 在 0 到 20 之间均等分成 200 份
t = linspace(0,20,200);
% 开环传递函数的分子多项式
num = 1;
for j=1:6
    % 求取开环传函的分母表达式
    den=conv([1,0], [T, 1+tau(j)]);
    % 建立开环传函
    s1 = tf(num*K, den);
    % 建立单位负反馈系统的传递函数
    sys = feedback(s1, 1);
    % 求取单位阶跃响应
    y(:,j) = step(sys, t);
end
% 在同一图上绘制单位阶跃响应曲线并添加栅格
plot(t, y(:,1:6)); grid
% 放置 K 取不同值的文字注释
gtext('tau=0');
gtext('tau=0.05');
gtext('tau=0.2');
gtext('tau=0.5');
gtext('tau=1.0');
gtext('tau=2.4');

% Page119：已知闭环传函，求未知量取不同值的单位阶跃响应
clear;
clc;

% 3 个不同的微分时间常数
tau = [0, 0.2, 0.4];
% 将 t 在 0 到 8 之间均等分为 80 份
t = linspace(0, 8, 80);
% 开环传递函数的分子、分母多项式
num = 4;
den = [1 2 4];
for j=1:3
    % 建立系统传递函数
    sys = tf(conv(num, [tau(j), 1]),den);
    % 求取单位阶跃响应
    y(:,j) = step(sys, t);
end
% 将 3 条响应曲线绘制在同一个图傻姑娘并添加栅格
plot(t,y(:, 1:3)); grid;
title('比例微分控制，不同微分时间下的系统阶跃响应');
gtext('tau=0');
gtext('tau=0.2');
gtext('tau=0.4');

% Page119：已知传递函数，分析主导极点，并比较新系统与原系统的单位阶跃响应
clear;
clc;

% 原系统的增益
K = 147.3;
% 仿真时间
t = 0:0.1:6;
num0 = K*[1, 1.5];
% 传递函数分母的系数
den00 = [1, 2, 5];
den01 = [1, 10, 26];
den02 = [1, 1.7];
% 建立原系统传递函数
sys0 = tf(num0, conv(den00, conv(den01, den02)));
% 求原系统的阶跃响应
y0 = step(sys0, t);
num1 = 5;
% 建立主导极点所构成的系统传递函数
sys1 = tf(num1, den00);
% 求取主导极点所构成的系统的阶跃响应
y1 = step(sys1, t);
% 绘制阶跃响应曲线并添加栅格
plot(t, y0, t, y1); grid;
% 添加图标题
title('阶跃响应对比');
% 放置区别两条曲线的文字注释
gtext('原系统的单位阶跃响应');
gtext('主导极点构成的系统的单位阶跃响应');

% Page131：已知传递函数，判断稳定性，绘制零极点图
clear;
clc;

% 开环零点、极点、增益
z = -2.5;
p = [0, -0.5, -0.7, -3];
k = 0.2;
% 建立零极点形式的开环传递函数
Go = zpk(z,p,k);
% 单位负反馈连接
Gc = feedback(Go, 1);
% 建立闭环传递函数
Gctf = tf(Gc);
% 获取闭环传递函数的特征多项式
dc = Gctf.den;
% 将特征多项式系数转换为字符形式的函数，便于查看
dens = poly2str(dc{1}, 's');

% 获取其多项式系数
den = [1, 4.2, 3.95, 1.25, 0.5];
% 求取特征根
p = roots(den);
% 绘制零极点图
pzmap(Gctf);
grid

% Page132：建立模型，分析稳定性，绘制阶跃响应曲线，计算稳态误差，分析性能
clear;
clc;

% 传递函数的分子、分母多项式系数
num1 = [20];
den1 = [1 2 0];
% 二阶系统的传递函数
sys1 = tf(num1, den1);
% 微分环节传递函数的分子、分母多项式系数
num2 = [0.1 0];
den2 = [0 1];
% 微分环节的传递函数
sys2 = tf(num2, den2);
% 内环反馈的传递函数
sys_inner = feedback(sys1, sys2);
% 外环反馈的传递函数
sys_outer = feedback(sys_inner, 1);

% 根据求得的系统传递函数，利用 roots 命令判断系统的稳定性
% 闭环系统传递函数分母多项式系数
den = [1 4 20];
% 求闭环系统特征多项式的根
roots(den)
% 利用 pzmap 命令绘制系统的零极点图
pzmap(sys_outer);
grid on;

% 计算系统的阶跃响应
% 闭环系统传递函数分子、分母多项式系数
num = [20];
den = [1 4 20];
% 计算闭环系统的阶跃响应
[y, t, x] = step(num, den);
% 绘制阶跃响应曲线
plot(x, y);
grid on;

% 计算系统的超调量
% 阶跃响应的稳态值
y_stable = 1;
% 闭环系统阶跃响应的最大值
max_response = max(y);
% 阶跃响应的超调量
sigma = (max_response - y_stable) / y_stable;

% 计算系统的上升时间
% 遍历响应曲线
for i=1:length(y)
    % 如果某个时刻系统的输出值大于稳态值
    if y(i)>y_stable
        % 循环中断
        break;
    end
end
% 计算此时对应的时间，就是阶跃响应的上升时间
tr = x(i)
% 计算系统的峰值时间
% 查找系统阶跃响应的最大值
[max_response, index] = max(y);
% 计算此时对应的时间，就是阶跃响应的峰值时间
tp = x(index)
% 计算系统的调整时间 ---> 取物插带为 2
% 遍历响应曲线
for i=1:length(y)
    % 如果当前响应值在误差带内
    if max(y(i:length(y))) <= 1.02 * y_stable
        if min(y(i:length(y))) >= 0.98 * y_stable
            % 循环退出
            break;
        end
    end
end
% 计算此时对应的时间，就是系统阶跃响应的调整时间
ts = x(i)