画图函数

import numpy as np
import matplotlib.pyplot as plt
def plot_bode(magnitude, amplitude, phase):amplitude = 20*np.log10(amplitude)fig = plt.figure()fig.subplots_adjust(hspace=0.1, top=0.99, bottom=0.05, left=0.06, right=0.99)ax1 = fig.add_subplot(211)ax1.grid(True)ax1.plot(magnitude, amplitude)ax2 = fig.add_subplot(212)ax2.plot(magnitude, phase)ax2.set_yticks(np.linspace(180, -360, 13))ax2.set_ylim(phase.min()-10, phase.max()+10)ax2.grid(True)plt.show()
if __name__ == '__main__':index = np.linspace(-2, 3, 41)omega = np.power(10, index)phi = -(np.pi/2+np.arctan2(0.1*omega, 1)+np.arctan2(0.2*omega, 1))*180/np.piamp = 30/(omega*np.sqrt(1+(0.1*omega)**2)*np.sqrt(1+(0.2*omega)**2))amp = 20*np.log10(amp)plot_bode(index, amp, phi)

一阶系统

G(s)=Ks+aG(s)=\frac{K}{s+a}G(s)=s+aK

import numpy as np
from zhnbodeplot import plot_bode
a = 10
K = 10
index = np.linspace(-1, 3, 41)
omega = np.power(10, index)
phi = -(np.pi/2+np.arctan2(0.1*omega, 1)+np.arctan2(0.2*omega, 1))*180/np.pi
amp = 30/(omega*np.sqrt(1+(0.1*omega)**2)*np.sqrt(1+(0.2*omega)**2))
plot_bode(index, amp, phi)

二阶系统

G(s)=ωn2s2+2ξωns+ωn2G(s)=\frac{\omega_n^2}{s^2+2\xi\omega_ns+\omega_n^2}G(s)=s2+2ξωns+ωn2ωn2

import numpy as np
from zhnbodeplot import plot_bode
omegan = 10
xi = 0.5
index = np.linspace(-1, 3, 41)
omega = np.power(10, index)
amp = omegan**2/np.sqrt((omegan**2 - omega**2)**2 + (2*xi*omegan*omega)**2)
amp = 20*np.log10(amp)
phi = -np.arctan2(2*xi*omegan*omega, omegan**2 - omega**2)
phi *= 180/np.pi
plot_bode(index, amp, phi)

时滞系统

G(s)=2e−τss+1G(s)=\frac{2\text{e}^{-\tau s}}{s+1}G(s)=s+12e−τs

import numpy as np
from zhnbodeplot import plot_bode
tau = 2*np.pi/3/np.sqrt(3)
index = np.linspace(-1, 0.5, 41)
omega = np.power(10, index)
amp = 2/np.sqrt(omega**2 + 1)
phi = -omega*tau-np.arctan2(omega, 1)
amp = 20*np.log10(amp)
phi *= 180/np.pi
plot_bode(index, amp, phi)

例

G(s)=30s(1+0.1s)(1+0.2s)G(s)=\frac{30}{s(1+0.1s)(1+0.2s)}G(s)=s(1+0.1s)(1+0.2s)30

phi = -(np.pi/2+np.arctan2(0.1*omega, 1)+np.arctan2(0.2*omega, 1))*180/np.pi
amp = 30/(omega*np.sqrt(1+(0.1*omega)**2)*np.sqrt(1+(0.2*omega)**2))

PID

G(s)=Kp+Kis+KdsG(s)=K_p+\frac{K_i}{s}+K_dsG(s)=Kp+sKi+Kds

Kp = 5
Ki = 1
Kd = 0.01
Re1 = Kp
Im1 = Kd*omega - Ki/omega
phi = np.arctan2(Im1, Re1)*180/np.pi
amp = np.sqrt(Re1**2 + Im1**2)

带滤波的PID

G(s)=Kp+Kis+Kdss+N=(Kp+Kd)s2+(KpN+Ki)s+KiNs2+Ns=KiN−(Kp+Kd)ω2+jω(KpN+Ki)−ω2+jNωG(s)=K_p+\frac{K_i}{s}+\frac{K_ds}{s+N}=\frac{(K_p+K_d)s^2+(K_pN+K_i)s+K_iN}{s^2+Ns} \\ =\frac{K_iN-(K_p+K_d)\omega^2+j\omega(K_pN+K_i)}{-\omega^2+jN\omega} G(s)=Kp+sKi+s+NKds=s2+Ns(Kp+Kd)s2+(KpN+Ki)s+KiN=−ω2+jNωKiN−(Kp+Kd)ω2+jω(KpN+Ki)

Kp = 5
Ki = 1
Kd = 0.01
N = 100
Re1 = Ki*N-(Kp+Kd)*omega**2
Im1 = omega*(Kp*N+Ki)
Re2 = -omega**2
Im2 = N*omega
phi = np.arctan2(Im1, Re1)*180/np.pi
phi -= np.arctan2(Im2, Re2)*180/np.pi
amp = np.sqrt(Re1**2 + Im1**2)
amp /= np.sqrt(Re2**2 + Im2**2)

二阶系统的谐振频率和峰值

A=ωn2(ωn2−ω2)2+(2ξωnω)2A=\frac{\omega_n^2}{\sqrt{(\omega_n^2-\omega^2)^2+(2\xi\omega_n\omega)^2}}A=(ωn2−ω2)2+(2ξωnω)2ωn2
令x=ω2x=\omega^2x=ω2，a=ωn2a=\omega_n^2a=ωn2，k+2=4ξ2k+2=4\xi^2k+2=4ξ2，
f(x)=x2−2ax+a2+(k+2)axf(x)=x^2-2ax+a^2+(k+2)axf(x)=x2−2ax+a2+(k+2)ax
f(x)f(x)f(x)最小值
f(−ka2)=a2−k2a24f(-\frac{ka}{2})=a^2-\frac{k^2a^2}{4}f(−2ka)=a2−4k2a2
x>0x>0x>0，a>0a>0a>0，所以k<0k<0k<0，即4ξ2−2<04\xi^2-2<04ξ2−2<0
ωr=−ωn2(4ξ2−2)2=ωn1−2ξ2\omega_r=\sqrt{-\frac{\omega_n^2(4\xi^2-2)}{2}}=\omega_n\sqrt{1-2\xi^2}ωr=−2ωn2(4ξ2−2)=ωn1−2ξ2
1Mr=1ωn2ωn4(1−(4ξ2−2)24)=1−(2ξ−1)2=2ξ1−ξ2\frac{1}{M_r}=\frac{1}{\omega_n^2}\sqrt{\omega_n^4(1-\frac{(4\xi^2-2)^2}{4})}=\sqrt{1-(2\xi-1)^2}=2\xi\sqrt{1-\xi^2}Mr1=ωn21ωn4(1−4(4ξ2−2)2)=1−(2ξ−1)2=2ξ1−ξ2
4x2(1−x2)≤14x^2(1-x^2)\le 14x2(1−x2)≤1，所以谐振峰值总是≥1。

二阶系统的幅相曲线

G(s)=ωn2s(s+2ξωn)G(s)=\frac{\omega_n^2}{s(s+2\xi\omega_n)}G(s)=s(s+2ξωn)ωn2

import matplotlib.pyplot as plt
import numpy as np
fig = plt.figure()
ax = plt.axes(projection='polar')index = np.linspace(-1, 3, 41)
radius2 = np.power(10, index)
theta2 =  np.pi/2 * np.ones(radius2.shape)
theta1 = np.linspace(-np.pi/2, np.pi/2, 100)
radius1 =  radius2[0] * np.ones(theta1.shape)
theta3 = np.linspace(np.pi/2, -np.pi/2, 100)
radius3 =  radius2[-1] * np.ones(theta1.shape)
stheta = np.concatenate([theta1, theta2, theta3, -theta2], axis=0)
sradius = np.concatenate([radius1, radius2, radius3, radius2[::-1]], axis=0)
# stheta = np.concatenate([-theta2], axis=0)
# sradius = np.concatenate([radius2], axis=0)xi = 0.5
omegan = 10
x = sradius * np.cos(stheta) + 2*xi*omegan
y = sradius * np.sin(stheta)
gradius = sradius * np.sqrt(x**2 + y**2)
gradius = omegan**2 / gradius
gtheta = -stheta - np.arctan2(y, x)
# gradius = 20 * np.log10(gradius)
# mask = gradius < -70
# gradius = gradius[mask]
# gtheta = gtheta[mask]
ax.plot(gtheta, gradius)
plt.show()

整体分成4段，第一段是半径极小的半圆，第二段是虚轴上的0到正无穷，第三段是无穷大的半圆，第四段是虚轴上的负无穷到0。幅相曲线整体见下图：

第二段见下图（注意幅度，下图是整体图的局部，下同）：

第三段见下图

注意第三段有个不合常理的地方，那就是里面那个圆的顺序是逆时针，也就是说整个曲线并没有交叉，大概是下面这样：

如果ξ<0\xi<0ξ<0，整体图是下面这样：

python matplotlib绘制伯德图相关推荐

Python matplotlib 绘制量场图
复习回顾 matplotlib 是基于Python语言的开源项目,pyplot提供一系列绘制2D图形的方法.随着版本的迭代,matplotlib 模块也支持绘制3D图形mplot3d工具包,制作动态图 ...
不愧是摸鱼高手Python matplotlib 绘制频谱图都会，能怪老板不管
复习回顾 matplotlib 是Python专门用来绘制渲染的模块,其底层主要分为脚本层.美工层和后端.脚本层为我们提供常见图形绘制如折线.柱状.直方.饼图.以往文章这么详细的Python mat ...
超详细的Python matplotlib 绘制动态图
复习回顾在matplotlib模块中我们前面学习绘制如折线.柱状.散点.直方图等静态图形.我们都知道在matplotlib模块主要有三层脚本层为用户提供快捷的绘制图形方法,美工层接收到脚本层的命令后 ...
Python matplotlib 绘制散点图还不收藏起来
复习回顾我们在往前几期中对matplotlib模块学习,对常用的反映数据变化的折线图,对比数据类型差异的柱状图和反应数据频率分布情况的直方图. 往前内容快速查看超详细的Python matplot ...
Python matplotlib 绘制等高线图
前言我们在往期对matplotlib.pyplot()方法学习,到现在我们已经会绘制折线图.柱状图.散点等常规的图表啦(往期的内容如下,大家可以方便查看往期内容) python入门到进阶,爬虫数据分 ...
这么详细的Python matplotlib 绘制图形还不赶紧收藏
前言我们前面对matplotlib模块底层结构学习,对其pyplot类(脚本层)类提供的绘制折线图.柱状图.饼图.直方图等统计图表的相关方法,列举往期文章如下. 超详细的Python matplot ...
超详细的Python matplotlib 绘制柱状图
复习回顾 Python 为数据展示提供了大量优秀的功能包,其中 matplotlib 模块可以方便绘制制作折线图.柱状图.散点图等高质量的数据包. 关于 matplotlib 模块,我们前期已经对ma ...
Python matplotlib绘制雷达图
Python matplotlib绘制雷达图本篇文章介绍使用matplotlib绘制雷达图. 雷达图也被称为网络图,蜘蛛图,星图,蜘蛛网图,是一个不规则的多边形.雷达图可以形象地展示相同事物的多维指 ...
Python matplotlib绘制饼图
Python matplotlib绘制饼图继前面使用matplotlib绘制折线图.散点图.柱状图和直方图,本篇文章继续介绍使用matplotlib绘制饼图. 一.matplotlib绘制饼图 # ...

python matplotlib绘制伯德图

文章目录