Python 数学应用（三）（3）

Python 数学应用（三）（2）https://developer.aliyun.com/article/1506401

对平面图形进行三角剖分

正如我们在第三章中看到的，微积分和微分方程，我们经常需要将连续区域分解为更小、更简单的区域。在之前的示例中，我们将实数区间缩小为一系列长度较小的小区间。这个过程通常称为离散化。在本章中，我们正在处理二维图形，因此我们需要这个过程的二维版本。为此，我们将一个二维图形（在这个示例中是一个多边形）分解为一系列更小和更简单的多边形。所有多边形中最简单的是三角形，因此这是二维离散化的一个很好的起点。找到一组"铺砌"几何图形的三角形的过程称为三角剖分。

在这个示例中，我们将学习如何使用 Shapely 包对多边形（带有孔）进行三角剖分。

准备工作

在这个示例中，我们需要将 NumPy 包导入为np，将 Matplotlib 包导入为mpl，并将pyplot模块导入为plt。

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

我们还需要从 Shapely 包中获取以下项目：

from shapely.geometry import Polygon
from shapely.ops import triangulate

如何做…

以下步骤向您展示了如何使用 Shapely 包对带有孔的多边形进行三角剖分：

1. 首先，我们需要创建一个代表我们希望进行三角剖分的图形的Polygon对象：

polygon = Polygon(
    [(2.0, 1.0), (2.0, 1.5), (-4.0, 1.5), (-4.0, 0.5), 
       (-3.0, -1.5), (0.0, -1.5), (1.0, -2.0), (1.0, -0.5), 
         (0.0, -1.0), (-0.5, -1.0), (-0.5, 1.0)],
    holes=[np.array([[-1.5, -0.5], [-1.5, 0.5], [-2.5, 0.5], 
       [-2.5, -0.5]])]
)

1. 现在，我们应该绘制图形，以便了解我们将在其中工作的区域：

fig, ax = plt.subplots()
plt_poly = mpl.patches.Polygon(polygon.exterior, 
   ec="k", lw="1", alpha=0.5, zorder=0)
ax.add_patch(plt_poly)
plt_hole = mpl.patches.Polygon(polygon.interiors[0], 
   ec="k", fc="w")
ax.add_patch(plt_hole)
ax.set(xlim=(-4.05, 2.05), ylim=(-2.05, 1.55))
ax.set_axis_off()

可以在下图中看到这个多边形。正如我们所看到的，这个图形中有一个"孔"，必须仔细考虑：

图 8.5：带有孔的示例多边形

1. 我们使用triangulate例程生成多边形的三角剖分。这个三角剖分包括外部边缘，这是我们在这个示例中不想要的：

triangles = triangulate(polygon)

1. 为了去除位于原始多边形外部的三角形，我们需要使用内置的filter例程，以及contains方法（在本章前面已经看到）：

filtered = filter(lambda p: polygon.contains(p), triangles) 

1. 将三角形绘制在原始多边形上，我们需要将 Shapely 三角形转换为 Matplotlib Patch对象，然后将其存储在PatchCollection中：

patches = map(lambda p: mpl.patches.Polygon(p.exterior), filtered)
col = mpl.collections.PatchCollection(patches, fc="none", ec="k")

1. 最后，我们将三角形补丁的集合添加到之前创建的图形中：

ax.add_collection(col)

在原始多边形上绘制的三角剖分可以在下图中看到。在这里，我们可以看到每个顶点都连接到另外两个顶点，形成了覆盖整个原始多边形的三角形系统：

图 8.6：带有孔的示例多边形的三角剖分

它是如何工作的…

triangulate例程使用一种称为Delaunay 三角剖分的技术将一组点连接到一组三角形中。在这种情况下，这组点是多边形的顶点。Delaunay 方法以这样一种方式找到这些三角形，即没有任何点包含在任何三角形的外接圆内。这是该方法的技术条件，但这意味着三角形被有效地选择，因为它避免了非常长、细的三角形。得到的三角剖分利用了原始多边形中存在的边缘，并连接了一些外部边缘。

为了去除原多边形外的三角形，我们使用内置的filter例程，它通过移除标准函数失败的项目来创建一个新的可迭代对象。这与 Shapely Polygon对象上的contains方法一起使用，以确定每个三角形是否位于原始图形内。正如我们之前提到的，我们需要将这些 Shapely 项目转换为 Matplotlib 补丁，然后才能将它们添加到图中。

还有更多…

三角剖分通常用于将复杂的几何图形简化为一组三角形，这些三角形对于某种计算任务来说要简单得多。然而，它们也有其他用途。三角剖分的一个特别有趣的应用是解决“艺术画廊问题”。这个问题涉及找到必要的“守卫”艺术画廊的最大数量。三角剖分是 Fisk 对艺术画廊定理的简单证明的重要部分，这个定理最初是由 Chvátal 证明的。

假设这个食谱中的多边形是一个艺术画廊的平面图，并且一些守卫需要放置在顶点上。一点工作就会表明，你需要在多边形的顶点处放置三个守卫，整个博物馆才能被覆盖。在下面的图像中，我们绘制了一个可能的布局：

图 8.7：在顶点上放置守卫的艺术画廊问题的一个可能解决方案。

点由点表示，并且它们相应的视野范围被阴影表示。

每个顶点都放置了一个守卫，并且他们的视野范围由相应的阴影区域表示。在这里，你可以看到整个多边形至少被一种颜色覆盖。艺术画廊问题的解决方案——实际上是原问题的一个变体——告诉我们，最多需要四名守卫。

另请参阅

关于艺术画廊问题的更多信息可以在 O’Rourke 的经典著作中找到：ORourke, J. (1987). Art gallery theorems and algorithms. New York: Oxford University Press.

计算凸包

如果图形内的每一对点都可以使用一条直线连接，并且这条直线也包含在图形内，那么几何图形被称为凸。凸体的简单例子包括点、直线、正方形、圆（圆盘）、正多边形等。图 8.5 中显示的几何图形不是凸的，因为孔的对面的点不能通过保持在图形内的直线连接起来。

从某种角度来看，凸图形是简单的，这意味着它们在各种应用中都很有用。一个特别的问题涉及找到包含一组点的最小凸集。这个最小凸集被称为这组点的凸包。

在这个食谱中，我们将学习如何使用 Shapely 包找到一组点的凸包。

准备工作

对于这个食谱，我们需要导入 NumPy 包作为np，导入 Matplotlib 包作为mpl，并导入plt模块：

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

我们还需要从 NumPy 导入默认的随机数生成器。我们可以这样导入：

from numpy.random import default_rng
rng = default_rng(12345)

最后，我们需要从 Shapely 导入MultiPoint类：

from shapely.geometry import MultiPoint

操作方法…

按照以下步骤找到一组随机生成点的凸包：

1. 首先，我们生成一个二维数组的随机数：

raw_points = rng.uniform(-1.0, 1.0, size=(50, 2))

1. 接下来，我们创建一个新图形，并在这个图形上绘制这些原始样本点：

fig, ax = plt.subplots()
ax.plot(raw_points[:, 0], raw_points[:, 1], "k.")
ax.set_axis_off()

这些随机生成的点可以在下图中看到。这些点大致分布在一个正方形区域内：

图 8.8：平面上的一组点

1. 接下来，我们构建一个MultiPoint对象，收集所有这些点并将它们放入一个单一对象中：

points = MultiPoint(raw_points)

1. 现在，我们使用convex_hull属性获取这个MultiPoint对象的凸包：

convex_hull = points.convex_hull

1. 然后，我们创建一个 MatplotlibPolygon补丁，可以在我们的图中绘制，以显示找到的凸包的结果：

patch = mpl.patches.Polygon(convex_hull.exterior, alpha=0.5,
   ec="k", lw=1.2)

1. 最后，我们将Polygon补丁添加到图中，以显示凸包：

ax.add_patch(patch)

随机生成的点的凸包可以在下图中看到：

图 8.9：平面上一组点的凸包

工作原理…

Shapely 包是围绕 GEOS 库的 Python 包装器，用于几何分析。Shapely 几何对象的convex_hull属性调用 GEOS 库中的凸包计算例程，从而产生一个新的 Shapely 对象。从这个教程中，我们可以看到一组点的凸包是一个多边形，其顶点是离“中心”最远的点。

构造贝塞尔曲线

贝塞尔曲线，或B 样条，是一族曲线，在矢量图形中非常有用-例如，它们通常用于高质量的字体包中。这是因为它们由少量点定义，然后可以用来廉价地计算沿曲线的大量点。这允许根据用户的需求来缩放细节。

在本教程中，我们将学习如何创建一个表示贝塞尔曲线的简单类，并计算沿其路径的若干点。

准备工作

在本教程中，我们将使用导入为np的 NumPy 包，导入为plt的 Matplotlib pyplot模块，以及 Python 标准库math模块中导入为binom的comb例程：

from math import comb as binom
import matplotlib.pyplot as plt
import numpy as np

如何做…

按照以下步骤定义一个表示贝塞尔曲线的类，该类可用于计算沿曲线的点：

1. 第一步是设置基本类。我们需要为实例属性提供控制点（节点）和一些相关的数字：

class Bezier:
    def __init__(self, *points):
        self.points = points
        self.nodes = n = len(points) - 1
        self.degree = l = points[0].size

1. 仍然在__init__方法中，我们生成贝塞尔曲线的系数，并将它们存储在实例属性的列表中：

self.coeffs = [binom(n, i)*p.reshape((l, 1)) for i, 
           p in enumerate(points)]

1. 接下来，我们定义一个__call__方法，使类可调用。我们将实例中的节点数加载到本地变量中，以便清晰明了：

def __call__(self, t):
        n = self.nodes

1. 接下来，我们重新整理输入数组，使其包含单行：

t = t.reshape((1, t.size))

1. 现在，我们使用实例的coeffs属性中的每个系数生成值数组的列表：

vals = [c @ (t**i)*(1-t)**(n-i) for i, 
           c in enumerate(self.coeffs)]

1. 最后，我们对步骤 5中构造的所有数组进行求和，并返回结果数组：

return np.sum(vals, axis=0)

1. 现在，我们将通过一个示例来测试我们的类。我们将为此示例定义四个控制点：

p1 = np.array([0.0, 0.0])
p2 = np.array([0.0, 1.0])
p3 = np.array([1.0, 1.0])
p4 = np.array([1.0, 3.0])

1. 接下来，我们为绘图设置一个新的图形，并用虚线连接线绘制控制点：

fig, ax = plt.subplots()
ax.plot([0.0, 0.0, 1.0, 1.0], [0.0, 1.0, 1.0, 3.0], "*--k")
ax.set(xlabel="x", ylabel="y", title="Bezier curve with 
    4 nodes, degree 3")

1. 然后，我们使用步骤 7中定义的四个点创建我们的Bezier类的新实例：

b_curve = Bezier(p1, p2, p3, p4)

1. 现在，我们可以使用linspace创建 0 到 1 之间等间距点的数组，并计算沿着贝塞尔曲线的点：

t = np.linspace(0, 1)
v = b_curve(t)

1. 最后，我们在之前绘制的控制点上绘制这条曲线：

ax.plot(v[0,:], v[1, :])

我们绘制的贝塞尔曲线可以在下图中看到。正如你所看到的，曲线从第一个点（0, 0）开始，结束于最终点（1, 3）：

图 8.10：使用四个节点构造的三次贝塞尔曲线

工作原理…

贝塞尔曲线由一系列控制点描述，我们以递归方式构造曲线。一个点的贝塞尔曲线是一个保持在该点的常数曲线。具有两个控制点的贝塞尔曲线是这两个点之间的线段：

当我们添加第三个控制点时，我们取对应点之间的线段，这些点是由一个较少点构成的贝塞尔曲线的曲线。这意味着我们使用以下公式构造具有三个控制点的贝塞尔曲线：

这种构造可以在下图中看到：

图 8.11：使用递归定义构造二次贝塞尔曲线。黑色虚线显示了两条线性贝塞尔曲线。

这种构造方式继续定义了任意数量控制点上的贝塞尔曲线。幸运的是，在实践中我们不需要使用这种递归定义，因为我们可以将公式展开成曲线的单一公式，即以下公式：

在这里，p[i]元素是控制点，t是一个参数，而

是二项式系数。请记住，t参数是生成曲线点的变化量。我们可以分离前述求和中涉及t的项和不涉及t的项。这定义了我们在步骤 2中定义的系数，每个系数由以下代码片段给出：

binom(n, i)*p.reshape((l, 1))

我们在这一步中对每个点p进行了 reshape，以确保它被排列为列向量。这意味着每个系数都是一个列向量（作为 NumPy 数组），由二项式系数缩放的控制点组成。

现在，我们需要指定如何在不同的t值上评估贝塞尔曲线。这就是我们利用 NumPy 包中的高性能数组操作的地方。在形成系数时，我们将控制点 reshape 为列向量。在步骤 4中，我们将输入t值 reshape 为行向量。这意味着我们可以使用矩阵乘法运算符将每个系数乘以相应的（标量）值，具体取决于输入的t。这就是步骤 5中列表推导式中发生的情况。在下一行中，我们将l×1数组乘以1×N数组，得到一个l×N数组：

c @ (t**i)*(1-t)**(n-i)

我们为每个系数都得到一个这样的数组。然后，我们可以使用np.sum例程来对这些l×N数组中的每一个进行求和，以得到贝塞尔曲线上的值。在本示例中，输出数组的顶行包含曲线的x值，底行包含曲线的y值。在指定axis=0关键字参数时，我们必须小心确保sum例程对我们创建的列表进行求和，而不是对该列表包含的数组进行求和。

我们定义的类是使用贝塞尔曲线的控制点进行初始化的，然后用于生成系数。曲线值的实际计算是使用 NumPy 完成的，因此这种实现应该具有相对良好的性能。一旦创建了这个类的特定实例，它的功能就非常像一个函数，正如你所期望的那样。但是，这里没有进行类型检查，所以我们只能用 NumPy 数组作为参数来调用这个“函数”。

还有更多…

贝塞尔曲线是使用迭代构造定义的，其中具有n个点的曲线是使用连接由第一个和最后一个n-1点定义的曲线来定义的。使用这种构造跟踪每个控制点的系数将很快导致我们用来定义前述曲线的方程。这种构造还导致贝塞尔曲线的有趣和有用的几何特性。

正如我们在这个配方的介绍中提到的，贝塞尔曲线出现在许多涉及矢量图形的应用程序中，比如字体。它们也出现在许多常见的矢量图形软件包中。在这些软件包中，通常会看到二次贝塞尔曲线，它们由三个点的集合定义。然而，你也可以通过提供两个端点以及这些点上的梯度线来定义一个二次贝塞尔曲线。这在图形软件包中更常见。生成的贝塞尔曲线将沿着梯度线离开每个端点，并在这些点之间平滑地连接曲线。

我们在这里构建的实现对于小型应用程序来说性能相对较好，但对于涉及在大量t值上渲染具有大量控制点的曲线的应用程序来说是不够的。对于这一点，最好使用一个用编译语言编写的低级软件包。例如，bezier Python 软件包使用编译的 Fortran 后端进行计算，并提供比我们在这里定义的类更丰富的接口。

当然，贝塞尔曲线可以自然地扩展到更高的维度。结果是一个贝塞尔曲面，使它们成为非常有用的通用工具，用于高质量、可伸缩的图形。

进一步阅读

• 计算几何中一些常见算法的描述可以在以下书籍中找到：Press, W.H., Teukolsky, S.A., Vetterling, W.T., and Flannery, B.P., 2007. Numerical recipes: the art of scientific computing*. 3rd ed. Cambridge: Cambridge University Press*。
  
• 有关计算几何中一些问题和技术的更详细描述，请查阅以下书籍：O’Rourke, J., 1994. Computational geometry in C*. Cambridge: Cambridge University Press*。
  

第十章：寻找最优解

在本章中，我们将讨论寻找给定情况下最佳结果的各种方法。这被称为优化，通常涉及最小化或最大化目标函数。目标函数是一个接受多个参数作为参数并返回代表给定参数选择的成本或回报的单个标量值的函数。关于最小化和最大化函数的问题实际上是相互等价的，因此我们只会在本章讨论最小化目标函数。最小化函数f(x)等同于最大化函数*-f*(x)。在我们讨论第一个配方时将提供更多细节。

我们可以利用的算法来最小化给定函数取决于函数的性质。例如，包含一个或多个变量的简单线性函数与具有许多变量的非线性函数相比，可用的算法不同。线性函数的最小化属于线性规划范畴，这是一个发展完善的理论。对于非线性函数，我们通常利用函数的梯度（导数）来寻找最小点。我们将讨论几种不同类型函数的最小化方法。

寻找单变量函数的极小值和极大值特别简单，如果函数的导数已知，可以轻松完成。如果不知道导数，则适用于适当配方的方法。最小化非线性函数配方中的注释提供了一些额外细节。

我们还将提供一个非常简短的介绍博弈论。广义上讲，这是一个围绕决策制定的理论，并在经济学等学科中具有广泛的影响。特别是，我们将讨论如何在 Python 中将简单的双人游戏表示为对象，计算与某些选择相关的回报，并计算这些游戏的纳什均衡。

我们将首先看如何最小化包含一个或多个变量的线性和非线性函数。然后，我们将继续研究梯度下降方法和使用最小二乘法进行曲线拟合。最后，我们将通过分析双人游戏和纳什均衡来结束本章。

在本章中，我们将涵盖以下配方：

• 最小化简单线性函数
  
• 最小化非线性函数
  
• 使用梯度下降法进行优化
  
• 使用最小二乘法拟合数据的曲线
  
• 分析简单的双人游戏
  
• 计算纳什均衡
  

让我们开始吧！

技术要求

在本章中，我们将像往常一样需要 NumPy 包、SciPy 包和 Matplotlib 包。我们还将需要 Nashpy 包用于最后两个配方。这些包可以使用您喜欢的包管理器（如pip）进行安装：

python3.8 -m pip install numpy scipy matplotlib nashpy

本章的代码可以在 GitHub 存储库的Chapter 09文件夹中找到，网址为github.com/PacktPublishing/Applying-Math-with-Python/tree/master/Chapter%2009。

查看以下视频以查看代码的实际操作：bit.ly/2BjzwGo。

最小化简单线性函数

在优化中我们面临的最基本问题是找到函数取得最小值的参数。通常，这个问题受到参数可能值的一些限制的约束，这增加了问题的复杂性。显然，如果我们要最小化的函数也很复杂，那么这个问题的复杂性会进一步增加。因此，我们必须首先考虑线性函数，它们的形式如下：

为了解决这些问题，我们需要将约束转化为计算机可用的形式。在这种情况下，我们通常将它们转化为线性代数问题（矩阵和向量）。一旦完成了这一步，我们就可以使用 NumPy 和 SciPy 中的线性代数包中的工具来找到我们所寻求的参数。幸运的是，由于这类问题经常发生，SciPy 有处理这种转化和随后求解的例程。

在这个配方中，我们将使用 SciPy optimize模块的例程来解决以下受限线性最小化问题：

这将受到以下条件的约束：

准备工作

对于这个配方，我们需要在别名np下导入 NumPy 包，以plt的名称导入 Matplotlib pyplot模块，以及 SciPy optimize模块。我们还需要从mpl_toolkits.mplot3d导入Axes3D类，以使 3D 绘图可用：

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

如何做…

按照以下步骤使用 SciPy 解决受限线性最小化问题：

1. 将系统设置为 SciPy 可以识别的形式：

A = np.array([
    [2, 1],   # 2*x0 + x1 <= 6
    [-1, -1]  # -x0 - x1 <= -4
])
b = np.array([6, -4])
x0_bounds = (-3, 14) # -3 <= x0 <= 14
x1_bounds = (2, 12)  # 2 <= x1 <= 12
c = np.array([1, 5])

1. 接下来，我们需要定义一个评估线性函数在x值处的例程，这是一个向量（NumPy 数组）：

def func(x):
    return np.tensordot(c, x, axes=1)

1. 然后，我们创建一个新的图并添加一组3d轴，我们可以在上面绘制函数：

fig = plt.figure()
ax = fig.add_subplot(projection="3d")
ax.set(xlabel="x0", ylabel="x1", zlabel="func")
ax.set_title("Values in Feasible region")

1. 接下来，我们创建一个覆盖问题区域的值网格，并在该区域上绘制函数的值：

X0 = np.linspace(*x0_bounds)
X1 = np.linspace(*x1_bounds)
x0, x1 = np.meshgrid(X0, X1)
z = func([x0, x1])
ax.plot_surface(x0, x1, z, alpha=0.3)

1. 现在，我们在函数值平面上绘制与临界线2*x0 + x1 == 6对应的线，并在我们的图上绘制落在范围内的值：

Y = (b[0] - A[0, 0]*X0) / A[0, 1]
I = np.logical_and(Y >= x1_bounds[0], Y <= x1_bounds[1])
ax.plot(X0[I], Y[I], func([X0[I], Y[I]]), "r", lw=1.5)

1. 我们重复这个绘图步骤，对第二条临界线x0 + x1 == -4：

Y = (b[1] - A[1, 0]*X0) / A[1, 1]
I = np.logical_and(Y >= x1_bounds[0], Y <= x1_bounds[1])
ax.plot(X0[I], Y[I], func([X0[I], Y[I]]), "r", lw=1.5)

1. 接下来，我们着色位于两条临界线之间的区域，这对应于最小化问题的可行区域：

B = np.tensordot(A, np.array([x0, x1]), axes=1)
II = np.logical_and(B[0, ...] <= b[0], B[1, ...] <= b[1]) 
ax.plot_trisurf(x0[II], x1[II], z[II], color="b", alpha=0.5)

函数值在可行区域上的图可以在以下图片中看到：

图 9.1：突出显示了可行区域的线性函数值

正如我们所看到的，位于这个着色区域内的最小值发生在两条临界线的交点处。

1. 接下来，我们使用linprog来解决带有我们在步骤 1中创建的边界的受限最小化问题。我们在终端中打印出结果对象：

res = optimize.linprog(c, A_ub=A, b_ub=b, bounds=
    (x0_bounds, x1_bounds))
print(res)

1. 最后，我们在可行区域上绘制最小函数值：

ax.plot([res.x[0]], [res.x[1]], [res.fun], "k*")

更新后的图可以在以下图片中看到：

图 9.2：在可行区域上绘制的最小值

在这里，我们可以看到linprog例程确实发现了最小值在两条临界线的交点处。

工作原理…

受限线性最小化问题在经济情况中很常见，您尝试在保持参数的其他方面的同时最小化成本。实际上，优化理论中的许多术语都反映了这一事实。解决这些问题的一个非常简单的算法称为单纯形法，它使用一系列数组操作来找到最小解。从几何上讲，这些操作代表着改变单纯形的不同顶点（我们在这里不定义），正是这一点赋予了算法其名称。

在我们继续之前，我们将简要概述单纯形法用于解决受限线性优化问题的过程。我们所面临的问题不是一个矩阵方程问题，而是一个矩阵不等式问题。我们可以通过引入松弛变量来解决这个问题，将不等式转化为等式。例如，通过引入松弛变量s[1]，可以将第一个约束不等式重写如下：

只要*s[1]不是负数，这就满足了所需的不等式。第二个约束不等式是大于或等于类型的不等式，我们必须首先将其更改为小于或等于类型。我们通过将所有项乘以-1 来实现这一点。这给了我们在配方中定义的矩阵A的第二行。引入第二个松弛变量s[2]*后，我们得到了第二个方程：

从中，我们可以构建一个矩阵，其列包含两个参数变量*x[1]和x[2]以及两个松弛变量s[1]和s[2]的系数。该矩阵的行代表两个边界方程和目标函数。现在可以使用对该矩阵进行初等行操作来解决这个方程组，以获得最小化目标函数的x[1]和x[2]*的值。由于解决矩阵方程很容易且快速，这意味着我们可以快速高效地最小化线性函数。

幸运的是，我们不需要记住如何将不等式系统化简为线性方程组，因为诸如linprog之类的例程会为我们完成这一工作。我们只需将边界不等式提供为矩阵和向量对，包括每个系数，以及定义目标函数的单独向量。linprog例程负责制定和解决最小化问题。

实际上，单纯形法并不是linprog例程用于最小化函数的算法。相反，linprog使用内点算法，这更有效率。（可以通过提供method关键字参数并使用适当的方法名称将方法设置为simplex或revised-simplex。在打印的输出结果中，我们可以看到只用了五次迭代就达到了解决方案。）该例程返回的结果对象包含最小值发生的参数值存储在x属性中，该最小值存储在fun属性中，以及有关解决过程的各种其他信息。如果方法失败，那么status属性将包含一个数值代码，描述方法失败的原因。

在本文的步骤 2中，我们创建了一个代表此问题的目标函数的函数。该函数以单个数组作为输入，该数组包含应在其上评估函数的参数空间值。在这里，我们使用了 NumPy 的tensordot例程（带有axes=1）来评估系数向量c与每个输入x的点积。在这里我们必须非常小心，因为我们传递给函数的值在后续步骤中将是一个 2×50×50 的数组。普通的矩阵乘法（np.dot）在这种情况下不会给出我们所期望的 50×50 数组输出。

在步骤 5和6中，我们计算了临界线上的点，这些点具有以下方程：

然后我们计算了相应的z值，以便绘制在由目标函数定义的平面上的线。我们还需要“修剪”这些值，以便只包括在问题中指定范围内的值。

还有更多…

本文介绍了受约束的最小化问题以及如何使用 SciPy 解决它。然而，相同的方法也可以用于解决受约束的最大化问题。这是因为最大化和最小化在某种意义上是对偶的，即最大化函数f(x)等同于最小化函数*-f*(x)，然后取其负值。事实上，我们在本文中使用了这一事实，将第二个约束不等式从≥改为≤。

在这个示例中，我们解决了一个只有两个参数变量的问题，但是相同的方法将适用于涉及两个以上这样的变量的问题（除了绘图步骤）。我们只需要向每个数组添加更多的行和列，以考虑这增加的变量数量 - 这包括提供给例程的边界元组。在处理非常大量的变量时，例程也可以与稀疏矩阵一起使用，以获得额外的效率。

linprog例程的名称来自线性规划，用于描述这种类型的问题 - 找到满足一些矩阵不等式的x的值，受其他条件的限制。由于与矩阵理论和线性代数的理论有非常紧密的联系，因此在线性规划问题中有许多非常快速和高效的技术，这些技术在非线性情境中是不可用的。

Python 数学应用（三）（4）https://developer.aliyun.com/article/1506406