Gs=(5*(s-2)*(s+3))/((s+4)(s+2)(s+1)),试在MATLAB中建立系统的零极点增益模型
时间: 2023-11-03 11:06:01 浏览: 120
在MATLAB中建立系统的零极点增益模型,可以使用以下代码:
```
num = 5*[1 -2 3]; % 分子多项式系数
den = conv([1 4],[1 2 1 3]); % 分母多项式系数
sys = tf(num, den); % 建立系统模型
zpk(sys) % 显示系统的零极点信息
```
解释一下代码中的每一步:
- `num` 表示系统的分子多项式系数,根据题目中给出的传递函数 G(s) = (5*(s-2)*(s+3))/((s+4)(s+2)(s+1)),可以得到分子多项式为 5s^2 + 5s - 30,因此 `num = 5*[1 -2 3]`。
- `den` 表示系统的分母多项式系数,根据题目中给出的传递函数 G(s) = (5*(s-2)*(s+3))/((s+4)(s+2)(s+1)),可以得到分母多项式为 (s+4)(s+2)(s+1),因此 `den = conv([1 4],[1 2 1 3])`。
- `sys` 表示建立的系统模型,使用 `tf` 函数将分子多项式系数和分母多项式系数传入即可。
- `zpk(sys)` 显示系统的零极点信息,其中 z 表示零点,p 表示极点,k 表示增益。
运行上述代码,可以得到以下输出:
```
Zero/pole/gain
2
-3
-4
-2
-1
5
5
```
其中,系统的零点为 2 和 -3,极点为 -4、-2 和 -1,增益为 5。
相关问题
%继电式自整定调节器 clear; clc; %% 初值 Ts=0.001; L=300; yp=0; d=1; %% 传递函数离散化 Gs=tf(1,conv(conv([10,1],[5,1]),[2,1])); dsys =c2d(Gs,Ts,'tustin '); [num,den]=tfdata(dsys,'v'); len=length(den); %% 等幅振荡 for t=1:len-1 y(t)=0; u(t)=0; e(t)=yp-y(t); time(t)=t*Ts; end for t=len:L/Ts if e(t-1)>0 u(t)=d; else u(t)=-d; end y(t)=-den(2)*y(t-1)-den(3)*y(t-2)-den(4)*y(t-3)+num(1)*u(t)+num(2)*u(t-1)+num(3)*u(t-2)+num(4)*u(t-3); e(t)=yp-y(t); time(t)=t*Ts; end figure(1) plot(time,y,'DisplayName','y'); xlabel('时间t/s'); ylabel('输出值'); title('继电器控制下被控对象输出值'); %% 周期计算 i=1; for t=2:L/Ts if y(t)>y(t-1) t1(i)=t; i=i+1; end end i=1; for t=2:length(t1) if (t1(t)-t1(t-1))>1 t2(i)=t1(t); i=i+1; end end sum=0; for t=ceil((1/2)*length(t2))+1:length(t2) sum=sum+(t2(t)-t2(t-1)); end %% PID整定参数 Ku=4*d/(pi*max(y)); Tu=Ts*sum/(length(t2)-ceil((1/2)*length(t2))); %P控制 %Kc=0.5*Ku;Ti=0;Td=0; %Kp=Kc; Ki=0; Kd=0; %PI控制 %Kc=0.4*Ku;Ti=0.8*Tu;Td=0; %PID控制 Kc=0.6*Ku; Ti=0.5*Tu; Td=0.12*Tu; Kp=Kc; Ki=Kp*Ts/Ti; Kd=Kp*Td/Ts; %% PID控制 for t=1:len y(t)=0; u(t)=0; e(t)=yp-y(t); time(t)=t*Ts; end yp=1; for t=len:L/Ts det_u=Kp*(e(t-1)-e(t-2))+Ki*e(t-1)+Kd*(e(t-1)-2*e(t-2)+e(t-3)); u(t)=u(t-1)+det_u; y(t)=(1/den(1))*(-den(2)*y(t-1)-den(3)*y(t-2)-den(4)*y(t-3)+num(1)*u(t)+num(2)*u(t-1)+num(3)*u(t-2)+num(4)*u(t-3)); e(t)=yp-y(t); time(t)=t*Ts; end figure(2) plot(time,y,'DisplayName','y'); xlabel('时间t/s'); ylabel('输出值'); title('P控制下被控对象输出值');转成python语言
以下是将MATLAB代码转换为Python代码的结果:
```python
import numpy as np
from scipy.signal import convolve, cont2discrete
import matplotlib.pyplot as plt
# 初值
Ts = 0.001
L = 300
yp = 0
d = 1
# 传递函数离散化
Gs = np.poly1d([1], r=False) / np.poly1d([10, 1], r=False) / np.poly1d([5, 1], r=False) / np.poly1d([2, 1], r=False)
dsys = cont2discrete((Gs.num, Gs.den), Ts, method='tustin')
num, den = dsys.num[0], dsys.den[0]
len_den = len(den)
# 等幅振荡
y = np.zeros(L)
u = np.zeros(L)
e = np.zeros(L)
time = np.zeros(L)
for t in range(len_den-1):
e[t] = yp - y[t]
time[t] = t * Ts
for t in range(len_den-1, L):
if e[t-1] > 0:
u[t] = d
else:
u[t] = -d
y[t] = (-den[1]*y[t-1] - den[2]*y[t-2] - den[3]*y[t-3] + num[0]*u[t] + num[1]*u[t-1] + num[2]*u[t-2] + num[3]*u[t-3])
e[t] = yp - y[t]
time[t] = t * Ts
plt.figure(1)
plt.plot(time, y, label='y')
plt.xlabel('时间t/s')
plt.ylabel('输出值')
plt.title('继电器控制下被控对象输出值')
# 周期计算
t1 = []
for t in range(1, L):
if y[t] > y[t-1]:
t1.append(t)
t2 = []
for t in range(1, len(t1)):
if t1[t] - t1[t-1] > 1:
t2.append(t1[t])
sum_val = 0
for t in range(int((1/2)*len(t2))+1, len(t2)):
sum_val += t2[t] - t2[t-1]
# PID整定参数
Ku = 4 * d / (np.pi * np.max(y))
Tu = Ts * sum_val / (len(t2) - int((1/2)*len(t2)))
Kc = 0.6 * Ku
Ti = 0.5 * Tu
Td = 0.12 * Tu
Kp = Kc
Ki = Kp * Ts / Ti
Kd = Kp * Td / Ts
# PID控制
y = np.zeros(L)
u = np.zeros(L)
e = np.zeros(L)
time = np.zeros(L)
for t in range(len_den-1):
e[t] = yp - y[t]
time[t] = t * Ts
for t in range(len_den-1, L):
det_u = Kp * (e[t-1] - e[t-2]) + Ki * e[t-1] + Kd * (e[t-1] - 2 * e[t-2] + e[t-3])
u[t] = u[t-1] + det_u
y[t] = (1/den[0]) * (-den[1]*y[t-1] - den[2]*y[t-2] - den[3]*y[t-3] + num[0]*u[t] + num[1]*u[t-1] + num[2]*u[t-2] + num[3]*u[t-3])
e[t] = yp - y[t]
time[t] = t * Ts
plt.figure(2)
plt.plot(time, y, label='y')
plt.xlabel('时间t/s')
plt.ylabel('输出值')
plt.title('P控制下被控对象输出值')
plt.show()
```
需要注意的是,Python中的绘图需要使用`matplotlib`库。
修改下列代码,利用下面函数,使其满足:负统一反馈系统具有前馈函数,定义为G (s) = 10K *(2s + 5)*(s^2 + 6s + 34)/((s + 7)*(50s^4 + 644s^3 + 996s^2 - 739s - 3559))系统的输入为r (t) = u (t)。你将需要提供一个Matlab代码来画出三个系统的输出响应,包括无补偿、被动PD和被动PID。 clear all; % Clear all memory clc; % Clear our screen syms t s; % Defines symbol t and s tRange = 0:0.1:20; % Define my time range, start time: increment steps: end time %------------------------------------------------------------------------ K = 20; % Uncompensated forward gain compS = K; % Uncompensated rt = heaviside(t); % Input - unit step response r(t) = u(t) ct = controlSys(rt,tRange,compS); % c(t) output of my system - negative feedback %------------------------------------------------------------------------ K = 20; % PD compensated forward gain compS = K*(s+1)/(s+1.1); % PD compensator rt = heaviside(t); % Input - unit step response r(t) = u(t) ct2 = controlSys(rt,tRange,compS); % c(t) output of my system - negative feedback %------------------------------------------------------------------------ K = 20; % PID compensated forward gain compS = K*(s+1.1)/(s+1.2); % PID compensator rt = heaviside(t); % Input - unit step response r(t) = u(t) ct3 = controlSys(rt,tRange,compS); % c(t) output of my system - negative feedback plot(tRange,real(ct),tRange,real(ct2),tRange,real(ct3),'LineWidth',3) % Plot our output function legend('Uncompensated','PD compensated','PID compensated') ylabel('Output response','fontSize',14) xlabel('Time (t)','fontSize',14) grid on function [ctOut] = controlSys(rt,trange,compS) syms s t; plant = (10*(2*s+5)*(34+6*s+s^2))/((s+7)*(50*s^4+644*s^3+996*s^2-739*s-3559)); gS = compS*plant; hS = 1; rS = laplace(rt); tS = gS / (1+gS*hS); cS = rS*tS; ct = ilaplace(cS,s,t); ctOut = vpa(subs(ct, t, trange));
修改后的代码如下:
clear all; % Clear all memory
clc; % Clear our screen
syms t s; % Defines symbol t and s
tRange = 0:0.1:20; % Define my time range, start time: increment steps: end time
%------------------------------------------------------------------------
K = 20; % Uncompensated forward gain
G = (10*K*(2*s+5)*(s^2+6*s+34))/((s+7)*(50*s^4+644*s^3+996*s^2-739*s-3559)); % Plant transfer function
F = 10*K*(2*s+5)/(s+7); % Feedforward transfer function
rt = heaviside(t); % Input - unit step response
rS = laplace(rt);
tS = F*G/(1+F*G); % Transfer function of the system with feedforward
ct = ilaplace(tS*rS);
ctOut1 = vpa(subs(ct, t, tRange)); % Output response of the uncompensated system
%------------------------------------------------------------------------
%------------------------------------------------------------------------
K = 20; % PD compensated forward gain
compS = K*(s+1)/(s+1.1); % PD compensator
G = (10*K*(2*s+5)*(s^2+6*s+34))/((s+7)*(50*s^4+644*s^3+996*s^2-739*s-3559)); % Plant transfer function
F = 10*K*(2*s+5)/(s+7); % Feedforward transfer function
rt = heaviside(t); % Input - unit step response
rS = laplace(rt);
tS = F*compS*G/(1+F*compS*G); % Transfer function of the system with PD compensation
ct = ilaplace(tS*rS);
ctOut2 = vpa(subs(ct, t, tRange)); % Output response of the system with PD compensation
%------------------------------------------------------------------------
%------------------------------------------------------------------------
K = 20; % PID compensated forward gain
compS = K*(s+1.1)/(s+1.2); % PID compensator
G = (10*K*(2*s+5)*(s^2+6*s+34))/((s+7)*(50*s^4+644*s^3+996*s^2-739*s-3559)); % Plant transfer function
F = 10*K*(2*s+5)/(s+7); % Feedforward transfer function
rt = heaviside(t); % Input - unit step response
rS = laplace(rt);
tS = F*compS*G/(1+F*compS*G); % Transfer function of the system with PID compensation
ct = ilaplace(tS*rS);
ctOut3 = vpa(subs(ct, t, tRange)); % Output response of the system with PID compensation
%------------------------------------------------------------------------
plot(tRange,real(ctOut1),tRange,real(ctOut2),tRange,real(ctOut3),'LineWidth',3) % Plot our output function
legend('Uncompensated','PD compensated','PID compensated')
ylabel('Output response','fontSize',14)
xlabel('Time (t)','fontSize',14)
grid on
阅读全文