手把手教你使用Numpy、Matplotlib、Scipy等5个Python库

导读：用Python设计控制系统，你会用到接下来要介绍的几个模块：Numpy、Matplotlib、Scipy、Sympy和Python-Control。

作者：南裕树（Yuki Minami）

来源：大数据DT（ID：hzdashuju）

01 Numpy

Numpy是用于数值计算的基础包。

使用Numpy可以高速高效地进行各种数值计算、统计处理以及信号处理。

import numpy as np

加载Numpy模块后，就可以使用平方根（np.sqrt）、绝对值（np.abs）、三角函数（np.sin、np.cos、np.arcsin、np.arccos）、指数（np.exp）、对数（np.log、np.log10）、四舍五入（np.round）等基本的数值计算用的函数了。

此外，弧度（rad）和角度（deg）换算的函数np.rad2deg也包括在内，还可以使用圆周率（np.pi）。

复数通过在虚部后面添加j来表示。可以使用np.imag来取出虚部。同样，可以使用np.real来取出实部，用np.conj来求得共轭复数。

还可以进行向量和矩阵运算。可以使用np.array来定义向量和矩阵。

A = np.array([ [1, 2], [-3, 4]])
print(A)

[[ 1  2]
 [-3  4]]

可以使用T方法来求得转置矩阵。

print(A.T)

[[ 1 -3]
 [ 2  4]]

可以使用np.linalg.inv来求得逆矩阵。

B = np.linalg.inv(A)
print(B)

[[ 0.4 -0.2]
 [ 0.3  0.1]]

可以像下面这样进一步求得矩阵对应的行列式的值、矩阵的秩、特征值以及范数。

np.linalg.det(A) #行列式

10.000000000000002

np.linalg.matrix_rank(A) #秩

2

w, v = np.linalg.eig(A) #特征值和特征向量
print('eigenvalue=',w)
print('eigenvector=\n',v)

eigenvalue= [2.5+1.93649167j 2.5-1.93649167j]
eigenvector=
 [[0.38729833-0.5j 0.38729833+0.5j]
 [0.77459667+0.j  0.77459667-0.j ]]

x = np.array([1, 2])
print(x)
np.linalg.norm(x)

2.23606797749979

此外，还可以使用np.arange构建数列。

Td = np.arange(0, 10, 1)
print(Td)

[0 1 2 3 4 5 6 7 8 9]

上例中构建了一个从0开始到10（不包含）为止，且间隔为1的数列。

02 Matplotlib

Matplotlib用于图形描绘，可以用它制作各种类型的图形和动画。绘图的具体例子可以参考它的主页:

http://matplotlib.org/gallery.html

可以像下面这样加载Matplotlib模块。

import matplotlib.pyplot as plt

让我们来看一个绘图的例子（见代码段2.1）。

• 代码段2.1 绘图

import numpy as np
import matplotlib.pyplot as plt
x = np.arange(0, 4 * np.pi, 0.1)
y = np.sin(x)
plt.plot(x, y) #横轴使用x，纵轴使用y绘图
plt.xlabel('x') #设定x轴的标签
plt.ylabel('y') #设定y轴的标签
plt.grid() #显示网格
plt.show()

这样可以画出如图2.9所示的图形。

▲图2.9 使用matplotlib绘图的例子

如果需要对细节进行进一步的调整以做出更美观的图形，则仅仅这样还是不够的。此时需要如代码段2.2所示采用能够对细节进行调整的面向对象的绘图方式。

• 代码段2.2 面向对象的绘图

fig, ax = plt.subplots() #生成Figure和Axes对象
ax.plot(x, y) #在Axes对象中生成图形
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.grid()
plt.show()

可以在Figure对象中生成Axes对象，并将数据传递给Axes对象的方法（.plot()）来进行绘图。可以进一步使用Axes对象的方法（.set_xlabel()等）来对图形进行调整。虽然可以像代码段2.1这样，在没有明确地生成对象的前提下描绘图形，但是如果明确生成各个对象，就可以在绘图的同时对细节进行调整。例如，可以如图2.10所示，在Figure中生成两个Axes对象。

▲图2.10 在Figure中生成两个Axes对象

可以使用代码段2.3生成图2.11。

• 代码段2.3 生成图2.11的代码

fig, ax = plt.subplots(2,1) #设置2行1列的子图形
 
x = np.arange(0, 4 * np.pi, 0.1)
y = np.sin(x)
z = np.cos(x)
w = y + z
 
#生成第一个图形
ax[0].plot(x, y, ls='-', label='sin', c='k')
ax[0].plot(x, z, ls='-.', label='cos', c='k')
ax[0].set_xlabel('x')
ax[0].set_ylabel('y, z')
ax[0].set_xlim(0, 4*np.pi)
ax[0].grid()
ax[0].legend()
 
#生成第二个图形
ax[1].plot(x, w, color='k', marker='.')
ax[1].set_xlabel('x')
ax[1].set_ylabel('w')
ax[1].set_xlim(0, 4*np.pi)
ax[1].grid(ls=':')
 
fig.tight_layout()

▲图2.11 绘制面向对象型的图形

图例可以通过在绘图时使用label = 'hogehoge'添加，并在最后通过ax.legend()输出。还可以使用类似x.legend(loc = 'best')这样的语句来指定图例的位置。指定位置时可以使用表2.1中的字符串或数值。

▼表2.1 指定图例位置时使用的字符串和数值

线的类型可以通过linestyle = '-'或者ls = '-'来设定。符号与线型的关系如下所示。

• -：实线　

• -.：点画线

• --：虚线

• ..：点线

线的粗细可以使用类似linewidth = 2或者lw = 2这样的代码来指定。线的颜色可以通过color = 'r'或者c = 'r'来指定。字符与颜色的关系如下所示。

• b：蓝

• g：绿

• r：红

• c：青

• m：品红

• y：黄

• k：黑

• w：白

记号的类型通过marker = 'o'来指定。部分记号的类型如下所示：

"."  ","  "o"  "v"  "^"  "<"  ">"  "1"  "2"  "3"  "4"  "8"

"s"  "p"  "*"  "h"  "H"  "+"  "x"  "D"  "d"  "|"  "_"  "$x$"

实际的显示效果（按照顺序）如图2.12所示。记号的大小可以像s = 10这样设定。

▲图2.12 记号的类型

可以使用fig.savefig("hogehoge.pdf")来保存图形。

03 Scipy

Scipy是用于数值计算算法的扩展包，提供了信号处理、优化和统计等功能的函数，还包含了控制系统分析和设计用的函数。这里介绍用于求解微分方程的odeint。

试求下述微分方程的数值积分：

假设输入为u(t)：

代码段2.4的输出结果如图2.13所示。通过def定义微分方程system，并将其与初始值y0和时间t一起传递给odeint。

• 代码段2.4 微分方程的数值积分

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
 
#定义微分方程
def system(y, t):
    if t < 10.0:
        u = 0.0
    else:
        u = 1.0
 
    dydt = (-y + u)/5.0 
    return dydt
 
#通过设定初始值和时间求解微分方程
y0 = 0.5
t = np.arange(0, 40, 0.04)
y = odeint(system, y0, t)
 
#绘图
fig, ax = plt.subplots()
ax.plot(t, y, label='y', c='k')
ax.plot(t, 1 * (t>=10), ls='--', label='u')
ax.set_xlabel('t')
ax.set_ylabel('y, u')
ax.legend(loc='best')
ax.grid(ls=':')

▲图2.13 微分方程的数值计算

04 Sympy

Sympy是用于符号计算的模块。通过它可以把变量当作符号处理并进行各种计算。例如，可以进行表达式展开、因式分解、微分、积分甚至是拉普拉斯变换。

import sympy as sp 
sp.init_printing()
s = sp.Symbol('s')
root = sp.solve(2 * s**2 +5*s+3, s)
print(root)

[-3/2, -1]

在上面的例子中，我们通过Symbol将变量定义为字符，使用solve求得方程的根。如果使用init_printing()，还可以输出LaTeX格式的结果。此时不要使用print，应直接输入root。

使用expand可以进行表达式展开，使用factor可以进行因式分解。

f = sp.expand( (s+1)*(s+2)**2, s)
print(f)

s**3 + 5*s**2 + 8*s + 4

g = sp.factor(f, s) 
print(g)

(s + 1)*(s + 2)**2

除了上面介绍的功能以外还有series（泰勒展开）、apart（部分分式分解）、laplace_transform（拉普拉斯变换）、inverse_laplace_transform（逆拉普拉斯变换）等功能。

05 Python-Control

本文通过实例介绍Python-Control函数。函数及其基本使用方法总结在表2.4～表2.8中，请读者在阅读过程中随时查阅。如果系统sys表示成以下形式：

那么num就表示[b_m,...,b₀]，den就表示[a_n,...a₀]，k就表示k，z就表示[z₁,...,z_m，p就表示[p₁,...,p_n]。如果系统sys表示成以下形式：

那么A就表示A，B就表示B，C就表示C，D就表示D。

可通过下述方法加载模块。

import control

也可以使用下述方法从模块中导入类MATLAB函数。

from control.matlab import *

▼表2.4 模型描述

▼表2.5 模型分析

▼表2.6 控制系统设计

▼表2.7 模型解析（from control import cannonical_form）

▼表2.8 控制系统设计（from control import mixsyn）

除了上面介绍的这些函数之外还有一些其他的方便实用的函数。

比如，求取根轨迹的函数rlocus、求取平衡降阶模型的函数balred、求解李雅普诺夫方程的函数lyap、用于H₂控制系统设计的函数h2syn以及用于H_∞控制系统设计的函数hinfsyn等。

详细内容可以参考Python Control Systems Library：

https://python-control.readthedocs.io/en/0.8.1/index.html

关于作者：南裕树（Yuki Minami），大阪大学大学院工学研究科机械工学专业教授，京都大学大学院情报学研究科博士。主要研究领域包含制造技术、系统工程、深度学习等。

本文摘编自《用Python轻松设计控制系统》，经出版方授权发布。（ISBN：9787111688112）

《用Python轻松设计控制系统》

点击上图了解及购买

转载请联系微信：DoctorData

推荐语：面向用Python进行控制系统设计的人，边学边用Python进行编程实践。重点阐述以状态空间模型为对象的现代控制系统的设计。全书分七章，包括： Pythhon环境的构建、数据和流控制、反馈控制、状态空间模型、闭环系统的控制设计、开环系统的控制设计、PID控制相位进延迟补偿，以及稳健性控制数字实现的高级控制系统设计。

刷刷视频