python 控制系统仿真库control （⼀）伯德图
1、安装
我使⽤的是pycharm，因此直接搜索control安装即可。注意，最好已经事先安装了scipy,numpy以及matplotlib这⼏个库以备不时之需。
2、⽂档地址
3、使⽤例⼦
1、创建控制系统：
import control as ctrl
sysTf=ctrl.tf([1],[1,2,1])
sysTf为由传递函数定义的系统。但是按照⽂档说明，这个库还⽀持⽤状态空间定义的系统。不过，作为⼀个⾃控的初学者，⽬前只⽤得到传递函数。
这个传递函数定义了⼀个振荡环节，相当于：
很多课本上都有这个环节的伯德图。若要画出这个图，则需要这么做：
import control as ctrl
import matplotlib.pyplot as plt
import numpy as np
psiList = [0.05,0.2,0.5,0.707,1.0]
sysTfList=[]取之不尽用之不竭
for psi in psiList:
sysTfList.append(ctrl.tf([1],[1,2*psi,1]))
mag,pha,omega=ctrl.bode(sysTfList,omega=np.logspace(start=-1,stop=1,num=200),dB=True,deg=Fal)
plt.show()
G (s )=
()+2s +1ωn s 2ωn ψ1
当然，图⽚还是有美中不⾜的地⽅：缺少图例。
4、适当修改源码中的⼀些部分
以下是我修改过的freqplot.py，在pycharm中按下ctrl+B追踪control.bode_plot()这个⽅法即可找到这个⽂件。
主要修改是增加了bode()的参数表，⽐原来的代码增加了⼀个legend变量，传⼊⼀个与sysList等长度的字符串列表作为图例的⽂字。原理是新建了⼀个名为lineList的列表，保存循环中sysList的每⼀项所绘出的线条对象。
⽬前仅修改了margins!=True情况时绘出的图形，因此在伯德图中可以显⽰图例。⽤法如下：
import control as ctrl
import matplotlib.pyplot as plt
import math
import numpy as np
psiList = [0.05,0.2,0.5,0.707,1.0]
sysTfList=[]
for psi in psiList:
sysTfList.append(ctrl.tf([1],[1,2*psi,1]))
三杰mag,pha,omega=ctrl.bode(sysTfList,omega=np.logspace(start=-1,stop=1,num=200),dB=True,deg=Fal,
娜读音legends=['0.05','0.2','0.5','0.707','1.0'])
关于母亲的古诗plt.show()
——————————————————————————————————————⽣成的图⽚如下：
以下是源代码：
# freqplot.py - frequency domain plots for control systems
#
# Author: Richard M. Murray
# Date: 24 May 09
#
# This file contains some standard control system plots: Bode plots,
# Nyquist plots and pole-zero diagrams. The code for Nichols charts
# is in nichols.py.
#
# Copyright (c) 2010 by California Institute of Technology
# All rights rerved.
#
# Redistribution and u in source and binary forms, with or without
# 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.
一分钟故事#
# 3. Neither the name of the California Institute of Technology nor
# the names of its contributors may be ud to endor or promote
# products derived from this software without specific prior
# written permission.
#
# 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 CALTECH
# OR THE 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.
#
# $Id$
import matplotlib
时间管理学
import matplotlib.pyplot as plt
import scipy as sp
import numpy as np
import math
from .ctrlutil import unwrap
from .bdalg import feedback
from .margins import stability_margins
__all__ = ['bode_plot', 'nyquist_plot', 'gangof4_plot',
'bode', 'nyquist', 'gangof4']
#
# Main plotting functions
#
# This ction of the code contains the functions for generating
# frequency domain plots
#
# Bode plot
def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
Plot=True, omega_limits=None, omega_num=None, margins=None,legends=[],*args, **kwargs): """
Bode plot for a system
Plots a Bode plot for the system over a (optional) frequency range.
Parameters
----------
子何恃而往翻译
syslist : linsys
List of linear input/output systems (single system is OK)
omega : list
List of frequencies in rad/c to be ud for frequency respon
dB : boolean
If True, plot result in dB
Hz : boolean
If True, plot frequency in Hz (omega must be provided in rad/c)
deg : boolean
If True, plot pha in degrees (el radians)
Plot : boolean
If True, plot magnitude and pha
omega_limits: tuple, list, ... of two values
Limits of the to generate frequency vector.
If Hz=True the limits are in Hz otherwi in rad/s.
omega_num: int
number of samples
margins : boolean
If True, plot gain and pha margin
\*args, \**kwargs:
Additional options to matplotlib (color, linestyle, etc)
Returns
-------
mag : array (list if len(syslist) > 1)
magnitude
pha : array (list if len(syslist) > 1)
pha : array (list if len(syslist) > 1)
pha in radians
omega : array (list if len(syslist) > 1)
frequency in rad/c
Notes
-----
1. Alternatively, you may u the lower-level method (mag, pha, freq)
= sys.freqresp(freq) to generate the frequency respon for a system,
but it returns a MIMO respon.
2. If a discrete time model is given, the frequency respon is plotted
along the upper branch of the unit circle, using the mapping z = exp(j
\omega dt) where omega ranges from 0 to pi/dt and dt is the discrete
timeba. If not timeba is specified (dt = True), dt is t to 1.
Examples
--------
>>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
>>> mag, pha, omega = bode(sys)
"""
# Set default values for options
from . import config
if dB is None:
dB = config.bode_dB
if deg is None:
deg = config.bode_deg
if Hz is None:
Hz = config.bode_Hz
# If argument was a singleton, turn it into a list
if not getattr(syslist, '__iter__', Fal):
syslist = (syslist,)
if omega is None:
if omega_limits is None:
# Select a default range if none is provided
omega = default_frequency_range(syslist, Hz=Hz, number_of_samples=omega_num) el:
omega_limits = np.array(omega_limits)
if Hz:
omega_limits *= 2. * math.pi
if omega_num:
omega = sp.logspace(np.log10(omega_limits[0]),
np.log10(omega_limits[1]),
num=omega_num,
endpoint=True)
el:
omega = sp.logspace(np.log10(omega_limits[0]),
np.log10(omega_limits[1]),
endpoint=True)
年终决算mags, phas, omegas, nyquistfrqs = [], [], [], []
lineList = []
for sys in syslist:
if sys.inputs > 1 or sys.outputs > 1:
# TODO: Add MIMO bode plots.
rai NotImplementedError("Bode is currently only implemented for SISO systems.") el:
omega_sys = np.array(omega)
if sys.isdtime(True):
nyquistfrq = 2. * math.pi * 1. / sys.dt / 2.
omega_sys = omega_sys[omega_sys < nyquistfrq]
# TODO: What distance to the Nyquist frequency is appropriate?
el:
nyquistfrq = None
