13518219792
def do_POST(self):
10年积累的网站设计、做网站经验，可以快速应对客户对网站的新想法和需求。提供各种问题对应的解决方案。让选择我们的客户得到更好、更有力的网络服务。我虽然不认识你，你也不认识我。但先网站设计后付款的网站建设流程，更有鄂伦春免费网站建设让你可以放心的选择与我们合作。
mpath,margs=urllib.splitquery(self.path)
datas = self.rfile.read(int(self.headers['content-length']))
self.do_action(mpath, datas)
def do_action(self, path, args):
self.outputtxt(path + args )
def outputtxt(self, content):
#指定返回编码
enc = "UTF-8"
content = content.encode(enc)
f = io.BytesIO()
f.write(content)
f.seek(0)
self.send_response(200)
self.send_header("Content-type", "text/html; charset=%s" % enc)
self.send_header("Content-Length", str(len(content)))
self.end_headers()
shutil.copyfileobj(f,self.wfile)
最优化
为什么要做最优化呢？因为在生活中，人们总是希望幸福值或其它达到一个极值，比如做生意时希望成本最小，收入最大，所以在很多商业情境中，都会遇到求极值的情况。
函数求根
这里「函数的根」也称「方程的根」，或「函数的零点」。
先把我们需要的包加载进来。import numpy as npimport scipy as spimport scipy.optimize as optimport matplotlib.pyplot as plt%matplotlib inline
函数求根和最优化的关系？什么时候函数是最小值或最大值？
两个问题一起回答：最优化就是求函数的最小值或最大值，同时也是极值，在求一个函数最小值或最大值时，它所在的位置肯定是导数为 0 的位置，所以要求一个函数的极值，必然要先求导，使其为 0，所以函数求根就是为了得到最大值最小值。
scipy.optimize 有什么方法可以求根？
可以用 scipy.optimize 中的 bisect 或 brentq 求根。f = lambda x: np.cos(x) - x # 定义一个匿名函数x = np.linspace(-5, 5, 1000) # 先生成 1000 个 xy = f(x) # 对应生成 1000 个 f(x)plt.plot(x, y); # 看一下这个函数长什么样子plt.axhline(0, color='k'); # 画一根横线，位置在 y=0
opt.bisect(f, -5, 5) # 求取函数的根0.7390851332155535plt.plot(x, y)plt.axhline(0, color='k')plt.scatter([_], [0], c='r', s=100); # 这里的 [_] 表示上一个 Cell 中的结果，这里是 x 轴上的位置，0 是 y 上的位置
求根有两种方法，除了上面介绍的 bisect，还有 brentq，后者比前者快很多。%timeit opt.bisect(f, -5, 5)%timeit opt.brentq(f, -5, 5)10000 loops, best of 3: 157 s per loopThe slowest run took 11.65 times longer than the fastest. This could mean that an intermediate result is being cached.10000 loops, best of 3: 35.9 s per loop
函数求最小化
求最小值就是一个最优化问题。求最大值时只需对函数做一个转换，比如加一个负号，或者取倒数，就可转成求最小值问题。所以两者是同一问题。
初始值对最优化的影响是什么？
举例来说，先定义个函数。f = lambda x: 1-np.sin(x)/xx = np.linspace(-20., 20., 1000)y = f(x)
当初始值为 3 值，使用 minimize 函数找到最小值。minimize 函数是在新版的 scipy 里，取代了以前的很多最优化函数，是个通用的接口，背后是很多方法在支撑。x0 = 3xmin = opt.minimize(f, x0).x # x0 是起始点，起始点最好离真正的最小值点不要太远plt.plot(x, y)plt.scatter(x0, f(x0), marker='o', s=300); # 起始点画出来，用圆圈表示plt.scatter(xmin, f(xmin), marker='v', s=300); # 最小值点画出来，用三角表示plt.xlim(-20, 20);
初始值为 3 时，成功找到最小值。
现在来看看初始值为 10 时，找到的最小值点。x0 = 10xmin = opt.minimize(f, x0).xplt.plot(x, y)plt.scatter(x0, f(x0), marker='o', s=300)plt.scatter(xmin, f(xmin), marker='v', s=300)plt.xlim(-20, 20);
由上图可见，当初始值为 10 时，函数找到的是局部最小值点，可见 minimize 的默认算法对起始点的依赖性。
那么怎么才能不管初始值在哪个位置，都能找到全局最小值点呢？
如何找到全局最优点？
可以使用 basinhopping 函数找到全局最优点，相关背后算法，可以看帮助文件，有提供论文的索引和出处。
我们设初始值为 10 看是否能找到全局最小值点。x0 = 10from scipy.optimize import basinhoppingxmin = basinhopping(f,x0,stepsize = 5).xplt.plot(x, y);plt.scatter(x0, f(x0), marker='o', s=300);plt.scatter(xmin, f(xmin), marker='v', s=300);plt.xlim(-20, 20);
当起始点在比较远的位置，依然成功找到了全局最小值点。
如何求多元函数最小值？
以二元函数为例，使用 minimize 求对应的最小值。def g(X): x,y = X return (x-1)**4 + 5 * (y-1)**2 - 2*x*yX_opt = opt.minimize(g, (8, 3)).x # (8,3) 是起始点print X_opt[ 1.88292611 1.37658521]fig, ax = plt.subplots(figsize=(6, 4)) # 定义画布和图形x_ = y_ = np.linspace(-1, 4, 100)X, Y = np.meshgrid(x_, y_)c = ax.contour(X, Y, g((X, Y)), 50) # 等高线图ax.plot(X_opt[0], X_opt[1], 'r*', markersize=15) # 最小点的位置是个元组ax.set_xlabel(r"$x_1$", fontsize=18)ax.set_ylabel(r"$x_2$", fontsize=18)plt.colorbar(c, ax=ax) # colorbar 表示颜色越深，高度越高fig.tight_layout()
画3D 图。from mpl_toolkits.mplot3d import Axes3Dfrom matplotlib import cmfig = plt.figure()ax = fig.gca(projection='3d')x_ = y_ = np.linspace(-1, 4, 100)X, Y = np.meshgrid(x_, y_)surf = ax.plot_surface(X, Y, g((X,Y)), rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)cset = ax.contour(X, Y, g((X,Y)), zdir='z',offset=-5, cmap=cm.coolwarm)fig.colorbar(surf, shrink=0.5, aspect=5);
曲线拟合
曲线拟合和最优化有什么关系？
曲线拟合的问题是，给定一组数据，它可能是沿着一条线散布的，这时要找到一条最优的曲线来拟合这些数据，也就是要找到最好的线来代表这些点，这里的最优是指这些点和线之间的距离是最小的，这就是为什么要用最优化问题来解决曲线拟合问题。
举例说明，给一些点，找到一条线，来拟合这些点。
先给定一些点：N = 50 # 点的个数m_true = 2 # 斜率b_true = -1 # 截距dy = 2.0 # 误差np.random.seed(0)xdata = 10 * np.random.random(N) # 50 个 x，服从均匀分布ydata = np.random.normal(b_true + m_true * xdata, dy) # dy 是标准差plt.errorbar(xdata, ydata, dy, fmt='.k', ecolor='lightgray');
上面的点整体上呈现一个线性关系，要找到一条斜线来代表这些点，这就是经典的一元线性回归。目标就是找到最好的线，使点和线的距离最短。要优化的函数是点和线之间的距离，使其最小。点是确定的，而线是可变的，线是由参数值，斜率和截距决定的，这里就是要通过优化距离找到最优的斜率和截距。
点和线的距离定义如下：def chi2(theta, x, y): return np.sum(((y - theta[0] - theta[1] * x)) ** 2)
上式就是误差平方和。
误差平方和是什么？有什么作用？
误差平方和公式为：
误差平方和大，表示真实的点和预测的线之间距离太远，说明拟合得不好，最好的线，应该是使误差平方和最小，即最优的拟合线，这里是条直线。
误差平方和就是要最小化的目标函数。
找到最优的函数，即斜率和截距。theta_guess = [0, 1] # 初始值theta_best = opt.minimize(chi2, theta_guess, args=(xdata, ydata)).xprint(theta_best)[-1.01442005 1.93854656]
上面两个输出即是预测的直线斜率和截距，我们是根据点来反推直线的斜率和截距，那么真实的斜率和截距是多少呢？-1 和 2，很接近了，差的一点是因为有噪音的引入。xfit = np.linspace(0, 10)yfit = theta_best[0] + theta_best[1] * xfitplt.errorbar(xdata, ydata, dy, fmt='.k', ecolor='lightgray');plt.plot(xfit, yfit, '-k');
最小二乘（Least Square）是什么？
上面用的是 minimize 方法，这个问题的目标函数是误差平方和，这就又有一个特定的解法，即最小二乘。
最小二乘的思想就是要使得观测点和估计点的距离的平方和达到最小，这里的“二乘”指的是用平方来度量观测点与估计点的远近（在古汉语中“平方”称为“二乘”），“最小”指的是参数的估计值要保证各个观测点与估计点的距离的平方和达到最小。
关于最小二乘估计的计算，涉及更多的数学知识，这里不想详述，其一般的过程是用目标函数对各参数求偏导数，并令其等于 0，得到一个线性方程组。具体推导过程可参考斯坦福机器学习讲义 第 7 页。def deviations(theta, x, y): return (y - theta[0] - theta[1] * x)theta_best, ier = opt.leastsq(deviations, theta_guess, args=(xdata, ydata))print(theta_best)[-1.01442016 1.93854659]
最小二乘 leastsq 的结果跟 minimize 结果一样。注意 leastsq 的第一个参数不再是误差平方和 chi2，而是误差本身 deviations，即没有平方，也没有和。yfit = theta_best[0] + theta_best[1] * xfitplt.errorbar(xdata, ydata, dy, fmt='.k', ecolor='lightgray');plt.plot(xfit, yfit, '-k');
非线性最小二乘
上面是给一些点，拟合一条直线，拟合一条曲线也是一样的。def f(x, beta0, beta1, beta2): # 首先定义一个非线性函数，有 3 个参数 return beta0 + beta1 * np.exp(-beta2 * x**2)beta = (0.25, 0.75, 0.5) # 先猜 3 个 betaxdata = np.linspace(0, 5, 50)y = f(xdata, *beta)ydata = y + 0.05 * np.random.randn(len(xdata)) # 给 y 加噪音def g(beta): return ydata - f(xdata, *beta) # 真实 y 和 预测值的差，求最优曲线时要用到beta_start = (1, 1, 1)beta_opt, beta_cov = opt.leastsq(g, beta_start)print beta_opt # 求到的 3 个最优的 beta 值[ 0.25525709 0.74270226 0.54966466]
拿估计的 beta_opt 值跟真实的 beta = (0.25, 0.75, 0.5) 值比较，差不多。fig, ax = plt.subplots()ax.scatter(xdata, ydata) # 画点ax.plot(xdata, y, 'r', lw=2) # 真实值的线ax.plot(xdata, f(xdata, *beta_opt), 'b', lw=2) # 拟合的线ax.set_xlim(0, 5)ax.set_xlabel(r"$x$", fontsize=18)ax.set_ylabel(r"$f(x, \beta)$", fontsize=18)fig.tight_layout()
除了使用最小二乘，还可以使用曲线拟合的方法，得到的结果是一样的。beta_opt, beta_cov = opt.curve_fit(f, xdata, ydata)print beta_opt[ 0.25525709 0.74270226 0.54966466]
有约束的最小化
有约束的最小化是指，要求函数最小化之外，还要满足约束条件，举例说明。
边界约束def f(X): x, y = X return (x-1)**2 + (y-1)**2 # 这是一个碗状的函数x_opt = opt.minimize(f, (0, 0), method='BFGS').x # 无约束最优化
假设有约束条件，x 和 y 要在一定的范围内，如 x 在 2 到 3 之间，y 在 0 和 2 之间。bnd_x1, bnd_x2 = (2, 3), (0, 2) # 对自变量的约束x_cons_opt = opt.minimize(f, np.array([0, 0]), method='L-BFGS-B', bounds=[bnd_x1, bnd_x2]).x # bounds 矩形约束fig, ax = plt.subplots(figsize=(6, 4))x_ = y_ = np.linspace(-1, 3, 100)X, Y = np.meshgrid(x_, y_)c = ax.contour(X, Y, f((X,Y)), 50)ax.plot(x_opt[0], x_opt[1], 'b*', markersize=15) # 没有约束下的最小值，蓝色五角星ax.plot(x_cons_opt[0], x_cons_opt[1], 'r*', markersize=15) # 有约束下的最小值，红色星星bound_rect = plt.Rectangle((bnd_x1[0], bnd_x2[0]), bnd_x1[1] - bnd_x1[0], bnd_x2[1] - bnd_x2[0], facecolor="grey")ax.add_patch(bound_rect)ax.set_xlabel(r"$x_1$", fontsize=18)ax.set_ylabel(r"$x_2$", fontsize=18)plt.colorbar(c, ax=ax)fig.tight_layout()
不等式约束
介绍下相关理论，先来看下存在等式约束的极值问题求法，比如下面的优化问题。
目标函数是 f(w)，下面是等式约束，通常解法是引入拉格朗日算子，这里使用 ββ 来表示算子，得到拉格朗日公式为
l 是等式约束的个数。
然后分别对 w 和ββ 求偏导，使得偏导数等于 0，然后解出 w 和βiβi，至于为什么引入拉格朗日算子可以求出极值，原因是 f(w) 的 dw 变化方向受其他不等式的约束，dw的变化方向与f(w)的梯度垂直时才能获得极值，而且在极值处，f(w) 的梯度与其他等式梯度的线性组合平行，因此他们之间存在线性关系。（参考《最优化与KKT条件》）
对于不等式约束的极值问题
常常利用拉格朗日对偶性将原始问题转换为对偶问题，通过解对偶问题而得到原始问题的解。该方法应用在许多统计学习方法中。有兴趣的可以参阅相关资料，这里不再赘述。def f(X): return (X[0] - 1)**2 + (X[1] - 1)**2def g(X): return X[1] - 1.75 - (X[0] - 0.75)**4x_opt = opt.minimize(f, (0, 0), method='BFGS').xconstraints = [dict(type='ineq', fun=g)] # 约束采用字典定义，约束方式为不等式约束，边界用 g 表示x_cons_opt = opt.minimize(f, (0, 0), method='SLSQP', constraints=constraints).xfig, ax = plt.subplots(figsize=(6, 4))x_ = y_ = np.linspace(-1, 3, 100)X, Y = np.meshgrid(x_, y_)c = ax.contour(X, Y, f((X, Y)), 50)ax.plot(x_opt[0], x_opt[1], 'b*', markersize=15) # 蓝色星星，没有约束下的最小值ax.plot(x_, 1.75 + (x_-0.75)**4, '', markersize=15)ax.fill_between(x_, 1.75 + (x_-0.75)**4, 3, color="grey")ax.plot(x_cons_opt[0], x_cons_opt[1], 'r*', markersize=15) # 在区域约束下的最小值ax.set_ylim(-1, 3)ax.set_xlabel(r"$x_0$", fontsize=18)ax.set_ylabel(r"$x_1$", fontsize=18)plt.colorbar(c, ax=ax)fig.tight_layout()
scipy.optimize.minimize 中包括了多种最优化算法，每种算法使用范围不同，详细参考官方文档。
多目标优化
目标优化问题一般地就是指通过一定的优化算法获得目标函数的最优化解。当优化的目标函数为一个时称之为单目标优化(Single-
objective Optimization Problem,
SOP)。当优化的目标函数有两个或两个以上时称为多目标优化(Multi-objective Optimization Problem,
MOP)。不同于单目标优化的解为有限解，多目标优化的解通常是一组均衡解。
多目标优化算法归结起来有传统优化算法和智能优化算法两大类。
1. 传统优化算法包括加权法、约束法和线性规划法等，实质上就是将多目标函数转化为单目标函数，通过采用单目标优化的方法达到对多目标函数的求解。
2. 智能优化算法包括进化算法（Evolutionary Algorithm, 简称EA）、粒子群算法（Particle Swarm Optimization, PSO）等。
Pareto最优解：
若x*∈C*,且在C中不存在比x更优越的解x，则称x*是多目标最优化模型式的Pareto最优解，又称为有效解。
一般来说，多目标优化问题并不存在一个最优解，所有可能的解都称为非劣解，也称为Pareto解。传统优化技术一般每次能得到Pareo解集中的一个，而
用智能算法来求解，可以得到更多的Pareto解，这些解构成了一个最优解集，称为Pareto最优解。它是由那些任一个目标函数值的提高都必须以牺牲其
他目标函数值为代价的解组成的集合,称为Pareto最优域，简称Pareto集。
Pareto有效(最优)解非劣解集是指由这样一些解组成的集合：与集合之外的任何解相比它们至少有一个目标函数比集合之外的解好。
求解多目标优化问题最有名的就是NSGA-II了，是多目标遗传算法，但其对解的选择过程可以用在其他优化算法上，例如粒子群，蜂群等等。这里简单介绍一下NSGA-II的选择算法。主要包含三个部分：
1. 快速非支配排序
要先讲一下支配的概念，对于解X1和X2，如果X1对应的所有目标函数都不比X2大（最小问题），且存在一个目标值比X2小，则X2被X1支配。
快速非支配排序是一个循环分级过程：首先找出群体中的非支配解集，记为第一非支配层，irank=1（irank是个体i的非支配值），将其从群体中除去，继续寻找群体中的非支配解集，然后irank=2。
2. 个体拥挤距离
为了使计算结果在目标空间比较均匀的分布，维持种群多样性，对每个个体计算拥挤距离，选择拥挤距离大的个体，拥挤距离的定义为：
L[i]d=L[i]d+(L[i+1]m−L[i−1]m)/(fmaxm−fminm)
L[i+1]m是第i+1个个体的第m目标函数值，fmaxm 和 fminm是集合中第m个目标函数的最大和最小值。
3. 精英策略选择
精英策略就是保留父代中的优良个体直接进入子代，防止获得的Pareto最优解丢失。将第t次产生的子代种群和父代种群合并，然后对合并后的新种群进行非支配排序，然后按照非支配顺序添加到规模为N的种群中作为新的父代。
想象一下，您有一个线性方程组和不等式系统。这样的系统通常有许多可能的解决方案。线性规划是一组数学和计算工具，可让您找到该系统的特定解，该解对应于某些其他线性函数的最大值或最小值。
混合整数线性规划是 线性规划 的扩展。它处理至少一个变量采用离散整数而不是连续值的问题。尽管乍一看混合整数问题与连续变量问题相似，但它们在灵活性和精度方面具有显着优势。
整数变量对于正确表示自然用整数表示的数量很重要，例如生产的飞机数量或服务的客户数量。
一种特别重要的整数变量是 二进制变量 。它只能取 零 或 一 的值，在做出是或否的决定时很有用，例如是否应该建造工厂或者是否应该打开或关闭机器。您还可以使用它们来模拟逻辑约束。
线性规划是一种基本的优化技术，已在科学和数学密集型领域使用了数十年。它精确、相对快速，适用于一系列实际应用。
混合整数线性规划允许您克服线性规划的许多限制。您可以使用分段线性函数近似非线性函数、使用半连续变量、模型逻辑约束等。它是一种计算密集型工具，但计算机硬件和软件的进步使其每天都更加适用。
通常，当人们试图制定和解决优化问题时，第一个问题是他们是否可以应用线性规划或混合整数线性规划。
以下文章说明了线性规划和混合整数线性规划的一些用例：
随着计算机能力的增强、算法的改进以及更多用户友好的软件解决方案的出现，线性规划，尤其是混合整数线性规划的重要性随着时间的推移而增加。
解决线性规划问题的基本方法称为，它有多种变体。另一种流行的方法是。
混合整数线性规划问题可以通过更复杂且计算量更大的方法来解决，例如，它在幕后使用线性规划。这种方法的一些变体是，它涉及使用 切割平面 ，以及。
有几种适用于线性规划和混合整数线性规划的合适且众所周知的 Python 工具。其中一些是开源的，而另一些是专有的。您是否需要免费或付费工具取决于问题的规模和复杂性，以及对速度和灵活性的需求。
值得一提的是，几乎所有广泛使用的线性规划和混合整数线性规划库都是以 Fortran 或 C 或 C++ 原生和编写的。这是因为线性规划需要对（通常很大）矩阵进行计算密集型工作。此类库称为求解器。Python 工具只是求解器的包装器。
Python 适合围绕本机库构建包装器，因为它可以很好地与 C/C++ 配合使用。对于本教程，您不需要任何 C/C++（或 Fortran），但如果您想了解有关此酷功能的更多信息，请查看以下资源：
基本上，当您定义和求解模型时，您使用 Python 函数或方法调用低级库，该库执行实际优化工作并将解决方案返回给您的 Python 对象。
几个免费的 Python 库专门用于与线性或混合整数线性规划求解器交互：
在本教程中，您将使用SciPy和PuLP来定义和解决线性规划问题。
在本节中，您将看到线性规划问题的两个示例：
您将在下一节中使用 Python 来解决这两个问题。
考虑以下线性规划问题：
你需要找到X和Ÿ使得红色，蓝色和黄色的不平等，以及不平等X 0和ÿ 0，是满意的。同时，您的解决方案必须对应于z的最大可能值。
您需要找到的自变量（在本例中为 x 和 y ）称为 决策变量 。要最大化或最小化的决策变量的函数（在本例中为 z） 称为 目标函数 、 成本函数 或仅称为 目标 。您需要满足的 不等式 称为 不等式约束 。您还可以在称为 等式约束 的约束中使用方程。
这是您如何可视化问题的方法：
红线代表的功能2 X + Ý = 20，和它上面的红色区域示出了红色不等式不满足。同样，蓝线是函数 4 x + 5 y = 10，蓝色区域被禁止，因为它违反了蓝色不等式。黄线是 x + 2 y = 2，其下方的黄色区域是黄色不等式无效的地方。
如果您忽略红色、蓝色和黄色区域，则仅保留灰色区域。灰色区域的每个点都满足所有约束，是问题的潜在解决方案。该区域称为 可行域 ，其点为 可行解 。在这种情况下，有无数可行的解决方案。
您想最大化z。对应于最大z的可行解是 最优解 。如果您尝试最小化目标函数，那么最佳解决方案将对应于其可行的最小值。
请注意，z是线性的。你可以把它想象成一个三维空间中的平面。这就是为什么最优解必须在可行区域的 顶点 或角上的原因。在这种情况下，最佳解决方案是红线和蓝线相交的点，稍后您将看到。
有时，可行区域的整个边缘，甚至整个区域，都可以对应相同的z值。在这种情况下，您有许多最佳解决方案。
您现在已准备好使用绿色显示的附加等式约束来扩展问题：
方程式 x + 5 y = 15，以绿色书写，是新的。这是一个等式约束。您可以通过向上一张图像添加相应的绿线来将其可视化：
现在的解决方案必须满足绿色等式，因此可行区域不再是整个灰色区域。它是绿线从与蓝线的交点到与红线的交点穿过灰色区域的部分。后一点是解决方案。
如果插入x的所有值都必须是整数的要求，那么就会得到一个混合整数线性规划问题，可行解的集合又会发生变化：
您不再有绿线，只有沿线的x值为整数的点。可行解是灰色背景上的绿点，此时最优解离红线最近。
这三个例子说明了 可行的线性规划问题 ，因为它们具有有界可行区域和有限解。
如果没有解，线性规划问题是 不可行的 。当没有解决方案可以同时满足所有约束时，通常会发生这种情况。
例如，考虑如果添加约束x + y 1会发生什么。那么至少有一个决策变量（x或y）必须是负数。这与给定的约束x 0 和y 0相冲突。这样的系统没有可行的解决方案，因此称为不可行的。
另一个示例是添加与绿线平行的第二个等式约束。这两行没有共同点，因此不会有满足这两个约束的解决方案。
一个线性规划问题是 无界的 ，如果它的可行区域是无界，将溶液不是有限。这意味着您的变量中至少有一个不受约束，可以达到正无穷大或负无穷大，从而使目标也无限大。
例如，假设您采用上面的初始问题并删除红色和黄色约束。从问题中删除约束称为 放松 问题。在这种情况下，x和y不会在正侧有界。您可以将它们增加到正无穷大，从而产生无限大的z值。
在前面的部分中，您研究了一个与任何实际应用程序无关的抽象线性规划问题。在本小节中，您将找到与制造业资源分配相关的更具体和实用的优化问题。
假设一家工厂生产四种不同的产品，第一种产品的日产量为x ₁，第二种产品的产量为x 2，依此类推。目标是确定每种产品的利润最大化日产量，同时牢记以下条件：
数学模型可以这样定义：
目标函数（利润）在条件 1 中定义。人力约束遵循条件 2。对原材料 A 和 B 的约束可以从条件 3 和条件 4 中通过对每种产品的原材料需求求和得出。
最后，产品数量不能为负，因此所有决策变量必须大于或等于零。
与前面的示例不同，您无法方便地将其可视化，因为它有四个决策变量。但是，无论问题的维度如何，原理都是相同的。
在本教程中，您将使用两个Python 包来解决上述线性规划问题：
SciPy 设置起来很简单。安装后，您将拥有开始所需的一切。它的子包 scipy.optimize 可用于线性和非线性优化。
PuLP 允许您选择求解器并以更自然的方式表述问题。PuLP 使用的默认求解器是COIN-OR Branch and Cut Solver (CBC)。它连接到用于线性松弛的COIN-OR 线性规划求解器 (CLP)和用于切割生成的COIN-OR 切割生成器库 (CGL)。
另一个伟大的开源求解器是GNU 线性规划工具包 (GLPK)。一些著名且非常强大的商业和专有解决方案是Gurobi、CPLEX和XPRESS。
除了在定义问题时提供灵活性和运行各种求解器的能力外，PuLP 使用起来不如 Pyomo 或 CVXOPT 等替代方案复杂，后者需要更多的时间和精力来掌握。
要学习本教程，您需要安装 SciPy 和 PuLP。下面的示例使用 SciPy 1.4.1 版和 PuLP 2.1 版。
您可以使用pip以下方法安装两者：
您可能需要运行pulptest或sudo pulptest启用 PuLP 的默认求解器，尤其是在您使用 Linux 或 Mac 时：
或者，您可以下载、安装和使用 GLPK。它是免费和开源的，适用于 Windows、MacOS 和 Linux。在本教程的后面部分，您将看到如何将 GLPK（除了 CBC）与 PuLP 一起使用。
在 Windows 上，您可以下载档案并运行安装文件。
在 MacOS 上，您可以使用 Homebrew：
在 Debian 和 Ubuntu 上，使用apt来安装glpk和glpk-utils：
在Fedora，使用dnf具有glpk-utils：
您可能还会发现conda对安装 GLPK 很有用：
安装完成后，可以查看GLPK的版本：
有关详细信息，请参阅 GLPK 关于使用Windows 可执行文件和Linux 软件包进行安装的教程。
在本节中，您将学习如何使用 SciPy优化和求根库进行线性规划。
要使用 SciPy 定义和解决优化问题，您需要导入scipy.optimize.linprog()：
现在您已经linprog()导入，您可以开始优化。
让我们首先解决上面的线性规划问题：
linprog()仅解决最小化（而非最大化）问题，并且不允许具有大于或等于符号 ( ) 的不等式约束。要解决这些问题，您需要在开始优化之前修改您的问题：
引入这些更改后，您将获得一个新系统：
该系统与原始系统等效，并且将具有相同的解决方案。应用这些更改的唯一原因是克服 SciPy 与问题表述相关的局限性。
下一步是定义输入值：
您将上述系统中的值放入适当的列表、元组或NumPy 数组中：
注意：请注意行和列的顺序！
约束左侧和右侧的行顺序必须相同。每一行代表一个约束。
来自目标函数和约束左侧的系数的顺序必须匹配。每列对应一个决策变量。
下一步是以与系数相同的顺序定义每个变量的界限。在这种情况下，它们都在零和正无穷大之间：
此语句是多余的，因为linprog()默认情况下采用这些边界（零到正无穷大）。
注：相反的float("inf")，你可以使用math.inf，numpy.inf或scipy.inf。
最后，是时候优化和解决您感兴趣的问题了。你可以这样做linprog()：
参数c是指来自目标函数的系数。A_ub和b_ub分别与不等式约束左边和右边的系数有关。同样，A_eq并b_eq参考等式约束。您可以使用bounds提供决策变量的下限和上限。
您可以使用该参数method来定义要使用的线性规划方法。有以下三种选择：
linprog() 返回具有以下属性的数据结构：
您可以分别访问这些值：
这就是您获得优化结果的方式。您还可以以图形方式显示它们：
如前所述，线性规划问题的最优解位于可行区域的顶点。在这种情况下，可行区域只是蓝线和红线之间的绿线部分。最优解是代表绿线和红线交点的绿色方块。
如果要排除相等（绿色）约束，只需删除参数A_eq并b_eq从linprog()调用中删除：
解决方案与前一种情况不同。你可以在图表上看到：
在这个例子中，最优解是红色和蓝色约束相交的可行（灰色）区域的紫色顶点。其他顶点，如黄色顶点，具有更高的目标函数值。
您可以使用 SciPy 来解决前面部分所述的资源分配问题：
和前面的例子一样，你需要从上面的问题中提取必要的向量和矩阵，将它们作为参数传递给.linprog()，然后得到结果：
结果告诉您最大利润是1900并且对应于x ₁ = 5 和x ₃ = 45。在给定条件下生产第二和第四个产品是没有利润的。您可以在这里得出几个有趣的结论：
opt.statusis0和opt.successis True，说明优化问题成功求解，最优可行解。
SciPy 的线性规划功能主要用于较小的问题。对于更大和更复杂的问题，您可能会发现其他库更适合，原因如下：
幸运的是，Python 生态系统为线性编程提供了几种替代解决方案，这些解决方案对于更大的问题非常有用。其中之一是 PuLP，您将在下一节中看到它的实际应用。
PuLP 具有比 SciPy 更方便的线性编程 API。您不必在数学上修改您的问题或使用向量和矩阵。一切都更干净，更不容易出错。
像往常一样，您首先导入您需要的内容：
现在您已经导入了 PuLP，您可以解决您的问题。
您现在将使用 PuLP 解决此系统：
第一步是初始化一个实例LpProblem来表示你的模型：
您可以使用该sense参数来选择是执行最小化（LpMinimize或1，这是默认值）还是最大化（LpMaximize或-1）。这个选择会影响你的问题的结果。
一旦有了模型，就可以将决策变量定义为LpVariable类的实例：
您需要提供下限，lowBound=0因为默认值为负无穷大。该参数upBound定义了上限，但您可以在此处省略它，因为它默认为正无穷大。
可选参数cat定义决策变量的类别。如果您使用的是连续变量，则可以使用默认值"Continuous"。
您可以使用变量x和y创建表示线性表达式和约束的其他 PuLP 对象：
当您将决策变量与标量相乘或构建多个决策变量的线性组合时，您会得到一个pulp.LpAffineExpression代表线性表达式的实例。
注意：您可以增加或减少变量或表达式，你可以乘他们常数，因为纸浆类实现一些Python的特殊方法，即模拟数字类型一样__add__()，__sub__()和__mul__()。这些方法用于像定制运营商的行为+，-和*。
类似地，您可以将线性表达式、变量和标量与运算符 ==、=以获取表示模型线性约束的纸浆.LpConstraint实例。
注：也有可能与丰富的比较方法来构建的约束.__eq__()，.__le__()以及.__ge__()定义了运营商的行为==，=。
考虑到这一点，下一步是创建约束和目标函数并将它们分配给您的模型。您不需要创建列表或矩阵。只需编写 Python 表达式并使用+=运算符将它们附加到模型中：
在上面的代码中，您定义了包含约束及其名称的元组。LpProblem允许您通过将约束指定为元组来向模型添加约束。第一个元素是一个LpConstraint实例。第二个元素是该约束的可读名称。
设置目标函数非常相似：
或者，您可以使用更短的符号：
现在您已经添加了目标函数并定义了模型。
注意：您可以使用运算符将 约束或目标附加到模型中，+=因为它的类LpProblem实现了特殊方法.__iadd__()，该方法用于指定 的行为+=。
对于较大的问题，lpSum()与列表或其他序列一起使用通常比重复+运算符更方便。例如，您可以使用以下语句将目标函数添加到模型中：
它产生与前一条语句相同的结果。
您现在可以看到此模型的完整定义：
模型的字符串表示包含所有相关数据：变量、约束、目标及其名称。
注意：字符串表示是通过定义特殊方法构建的.__repr__()。有关 的更多详细信息.__repr__()，请查看Pythonic OOP 字符串转换：__repr__vs__str__ .
最后，您已准备好解决问题。你可以通过调用.solve()你的模型对象来做到这一点。如果要使用默认求解器 (CBC)，则不需要传递任何参数：
.solve()调用底层求解器，修改model对象，并返回解决方案的整数状态，1如果找到了最优解。有关其余状态代码，请参阅LpStatus[]。
你可以得到优化结果作为 的属性model。该函数value()和相应的方法.value()返回属性的实际值：
model.objective持有目标函数model.constraints的值，包含松弛变量的值，以及对象x和y具有决策变量的最优值。model.variables()返回一个包含决策变量的列表：
如您所见，此列表包含使用 的构造函数创建的确切对象LpVariable。
结果与您使用 SciPy 获得的结果大致相同。
注意：注意这个方法.solve()——它会改变对象的状态，x并且y！
您可以通过调用查看使用了哪个求解器.solver：
输出通知您求解器是 CBC。您没有指定求解器，因此 PuLP 调用了默认求解器。
如果要运行不同的求解器，则可以将其指定为 的参数.solve()。例如，如果您想使用 GLPK 并且已经安装了它，那么您可以solver=GLPK(msg=False)在最后一行使用。请记住，您还需要导入它：
现在你已经导入了 GLPK，你可以在里面使用它.solve()：
该msg参数用于显示来自求解器的信息。msg=False禁用显示此信息。如果要包含信息，则只需省略msg或设置msg=True。
您的模型已定义并求解，因此您可以按照与前一种情况相同的方式检查结果：
使用 GLPK 得到的结果与使用 SciPy 和 CBC 得到的结果几乎相同。
一起来看看这次用的是哪个求解器：
正如您在上面用突出显示的语句定义的那样model.solve(solver=GLPK(msg=False))，求解器是 GLPK。
您还可以使用 PuLP 来解决混合整数线性规划问题。要定义整数或二进制变量，只需传递cat="Integer"或cat="Binary"到LpVariable。其他一切都保持不变：
在本例中，您有一个整数变量并获得与之前不同的结果：
Nowx是一个整数，如模型中所指定。（从技术上讲，它保存一个小数点后为零的浮点值。）这一事实改变了整个解决方案。让我们在图表上展示这一点：
如您所见，最佳解决方案是灰色背景上最右边的绿点。这是两者的最大价值的可行的解决方案x和y，给它的最大目标函数值。
GLPK 也能够解决此类问题。
现在你可以使用 PuLP 来解决上面的资源分配问题：
定义和解决问题的方法与前面的示例相同：
在这种情况下，您使用字典 x来存储所有决策变量。这种方法很方便，因为字典可以将决策变量的名称或索引存储为键，将相应的LpVariable对象存储为值。列表或元组的LpVariable实例可以是有用的。
上面的代码产生以下结果：
如您所见，该解决方案与使用 SciPy 获得的解决方案一致。最有利可图的解决方案是每天生产5.0第一件产品和45.0第三件产品。
让我们把这个问题变得更复杂和有趣。假设由于机器问题，工厂无法同时生产第一种和第三种产品。在这种情况下，最有利可图的解决方案是什么？
现在您有另一个逻辑约束：如果x ₁ 为正数，则x ₃ 必须为零，反之亦然。这是二元决策变量非常有用的地方。您将使用两个二元决策变量y ₁ 和y ₃，它们将表示是否生成了第一个或第三个产品：
除了突出显示的行之外，代码与前面的示例非常相似。以下是差异：
这是解决方案：
事实证明，最佳方法是排除第一种产品而只生产第三种产品。
就像有许多资源可以帮助您学习线性规划和混合整数线性规划一样，还有许多具有 Python 包装器的求解器可用。这是部分列表：
其中一些库，如 Gurobi，包括他们自己的 Python 包装器。其他人使用外部包装器。例如，您看到可以使用 PuLP 访问 CBC 和 GLPK。
您现在知道什么是线性规划以及如何使用 Python 解决线性规划问题。您还了解到 Python 线性编程库只是本机求解器的包装器。当求解器完成其工作时，包装器返回解决方案状态、决策变量值、松弛变量、目标函数等。
函数详见rres，此代码使该算法运行了两次
收获：
这是我第一个实现的代码。学习完该算法以后，逻辑框架基本上就有了，剩下需要明确的就是对应的python的语言。于是我就开始了查找“如何定义函数”（详见mofan的优酷），“循环体”和“if条件语句”的格式（）“数学符号”（详见mofan的优酷），以及print的使用
1.def是python中指定义，一般用来定义函数，如果需要深度学习搭建网络可用来定义网络。值得注意的一点是
我不清楚为什么，但是如果没有加的话，那个函数公式就是一个花瓶，就像一个结果输不出去。
2.最坑的就是逻辑。一开始逻辑没理清楚，或者说在代码上有疏漏，导致我将left和right放在了循环体里，结果可想而知。不过也是因为这个错误，我知道pycharm中的debug怎么用，挺简单的，百度一下就出来了。
3.不知道什么原因，看的莫烦视频中的print多个变量一起输出是没有办法在我的pycharm中使用的，出来的结果很奇怪。可能是因为我是win10不是ios吧。print如果多个变量一起输出必须是print("名字：%s,名字2：%s"%(a,b))结果输出就是名字：a ,名字2：b
关于python中数据变量。第一遍运行结果出现很明显不对，于是我采用了debug。结果发现，mid1处一直为1而不是1.5，于是就开始了解数据变量。起初我猜测python默认所有变量为整型，但是根据二分法的结果我意识到此猜测不对，所以要改整个file的变量格式没有必要。所以我就在mid1式子前面加了一个float，结果就显示为1.5了。但是如果我将整个式子用（）括起来，前面加float，结果还是1。我不太理解为什么。不过我知道了python的数据格式是根据输入量决定的，也就是说你的输入量如果是整型，那么与其直接相关的计算输出结果一定是整型，而且还是不采用进位的整型。在我没有采用+float/+.0这两种方法之前，mid1~3全部是整型。
或者不再mid1前面加float,直接将输入量后面点个点就行
真的很想吐槽一下print,好麻烦啊啊啊啊每次都得弄个%s,而且有时候还不能放一起！！！！
不要问我掌握了什么，要问我现在写完这个代码后有多么的爱python的精度表示 :-)我决定以后只要再编写数学公式的代码都将输入量的小数学点后面补很多0
fibonacci函数定义，每次debug后我的手都是抖的O( _ )O~
不知道自己什么时候有的强迫症，只要是代码下面有“~”我就必须要消掉。笑哭。这个很简单，前四个除了费波纳茨，都很简单。
这个公式看起来很麻烦，便写的时候更要谨慎。我上回把那个2搁在了分号下面，结果很大，所以还是换算成0.5更好（PS：勿忘那长河般的0）。
虽然代码很长，但是主要是因为print太多。本打算在开头print，最后结果会漏掉最后一部分。懒得想其他办法了，直接就这样吧
一开始while里面写成了,导致run不出来。继而，debug也没法用。在网上一查才知道 “没联网”+“没选断点”。最后想尝试将else里面的内容输出来，结果发现run以后被刷屏了。于是改成i7以后还是不行，于是想着加一个break跳出循环，结果成效了。
然后刚刚由debug了一下，才知道原来是i+1在if里面，因为没有办法+1，所以i=6一直存在，就不断循环。因为加break也好，i+1也好，都可以。
这是我第一组自己实现的python代码，就是数学公式用python语言组装起来。刚开始的时候知道大概需要在语言中体现什么，但不太清楚。于是我就在网上找了几个二分法的，他们都各有不同，但框架都差不多，不过如果要用到我们的那个公式里还需要改变很多。然后我就开始分析我们的题，我发现大体需要两部分，一部分函数定义，一部分循环体。但我不知道如何定义函数，如何写数学公式，如何弄变量，也就是说一些小点不太会，所以我选择直接百度。因为我知道自己阅读的能力不错，相比于从视频中提取要素，我更擅长通过阅读获得要点。有目的性地找知识点，掌握地更牢固。
于是我就开始了第一个——二分法的编写。我发现，自己出现了很多错误而且有很多地方都很基础。但我依然没选择视频，而是将这些问题直接在百度上找，因为视频讲完或许你也没找到点。当然，这是一步一步走的，不是直接就将程序摆上去，一点一点改。
随着前两个的成功，我发现自己对于这些代码有了自信，似乎看透了他们的伪装，抓住了本质。除此之外，我还意识到自己自从8月份以后，学习能力似乎提高了不少，而且有了更为有效的学习方法。各方面都有了一定的觉醒。除了第一个找了几个牛头不对马嘴的代码，其他都是根据自己的逻辑写，逻辑通下来以后，对应语言中某一部分不知道如何翻译就去百度，其实这几个套路都一样或者说数学公式转化的套路都一样。
我还意识到，汇编其实是最难的语言，目前为止所学到的，因为很多都需要自己去定义，去死抠，需要记住大量的指令且不能灵活变通。但是其他的却只需要将一些对应的记下来就好。python真的挺简单的。而且，我发现自己今天似乎打开了新世界的大门，我爱上了这种充满了灵性的东西，充满了严谨的美丽，还有那未知的变化，我发现我似乎爱上了代码。可能不仅仅局限于python，这些语言都充满了挑战性。我觉得当你疑惑的时候，就需要相信直觉，至少我发现它很准
0.导入需要的包import pandas as pd
import numpy as np
import statsmodels.api as sm #统计运算
import scipy.stats as scs #科学计算
import matplotlib.pyplot as plt #绘图
1.选取几只感兴趣的股票
000413 东旭光电，000063 中兴通讯，002007 华兰生物，000001 平安银行，000002 万科A
并比较一下数据（2015-01-01至2015-12-31）
In[1]:
stock_set = ['000413.XSHE','000063.XSHE','002007.XSHE','000001.XSHE','000002.XSHE']
noa = len(stock_set)
df = get_price(stock_set, start_date = '2015-01-01', end_date ='2015-12-31', 'daily', ['close'])
data = df['close']
#规范化后时序数据
(data/data.ix[0]*100).plot(figsize = (8,5))
Out[1]:
2.计算不同证券的均值、协方差
每年252个交易日，用每日收益得到年化收益。计算投资资产的协方差是构建资产组合过程的核心部分。运用pandas内置方法生产协方差矩阵。
In [2]:
returns = np.log(data / data.shift(1))
returns.mean()*252
Out[2]:
000413.XSHE 0.184516
000063.XSHE 0.176790
002007.XSHE 0.309077
000001.XSHE -0.102059
000002.XSHE 0.547441
In [3]:
returns.cov()*252
Out[3]:
3.给不同资产随机分配初始权重
由于A股不允许建立空头头寸，所有的权重系数均在0-1之间
In [4]:
weights = np.random.random(noa)
weights /= np.sum(weights)
weights
Out[4]:
array([ 0.37505798, 0.21652754, 0.31590981, 0.06087709, 0.03162758])
4.计算预期组合年化收益、组合方差和组合标准差
In [5]:
np.sum(returns.mean()*weights)*252
Out[5]:
0.21622558669017816
In [6]:
np.dot(weights.T, np.dot(returns.cov()*252,weights))
Out[6]:
0.23595133640121463
In [7]:
np.sqrt(np.dot(weights.T, np.dot(returns.cov()* 252,weights)))
Out[7]:
0.4857482232609962
5.用蒙特卡洛模拟产生大量随机组合
进行到此，我们最想知道的是给定的一个股票池（证券组合）如何找到风险和收益平衡的位置。
下面通过一次蒙特卡洛模拟，产生大量随机的权重向量，并记录随机组合的预期收益和方差。
In [8]:
port_returns = []
port_variance = []
for p in range(4000):
weights = np.random.random(noa)
weights /=np.sum(weights)
port_returns.append(np.sum(returns.mean()*252*weights))
port_variance.append(np.sqrt(np.dot(weights.T, np.dot(returns.cov()*252, weights))))
port_returns = np.array(port_returns)
port_variance = np.array(port_variance)
#无风险利率设定为4%
risk_free = 0.04
plt.figure(figsize = (8,4))
plt.scatter(port_variance, port_returns, c=(port_returns-risk_free)/port_variance, marker = 'o')
plt.grid(True)
plt.xlabel('excepted volatility')
plt.ylabel('expected return')
plt.colorbar(label = 'Sharpe ratio')
Out[8]:
6.投资组合优化1——sharpe最大
建立statistics函数来记录重要的投资组合统计数据（收益，方差和夏普比）
通过对约束最优问题的求解，得到最优解。其中约束是权重总和为1。
In [9]:
def statistics(weights):
weights = np.array(weights)
port_returns = np.sum(returns.mean()*weights)*252
port_variance = np.sqrt(np.dot(weights.T, np.dot(returns.cov()*252,weights)))
return np.array([port_returns, port_variance, port_returns/port_variance])
#最优化投资组合的推导是一个约束最优化问题
import scipy.optimize as sco
#最小化夏普指数的负值
def min_sharpe(weights):
return -statistics(weights)[2]
#约束是所有参数(权重)的总和为1。这可以用minimize函数的约定表达如下
cons = ({'type':'eq', 'fun':lambda x: np.sum(x)-1})
#我们还将参数值(权重)限制在0和1之间。这些值以多个元组组成的一个元组形式提供给最小化函数
bnds = tuple((0,1) for x in range(noa))
#优化函数调用中忽略的唯一输入是起始参数列表(对权重的初始猜测)。我们简单的使用平均分布。
opts = sco.minimize(min_sharpe, noa*[1./noa,], method = 'SLSQP', bounds = bnds, constraints = cons)
opts
Out[9]:
status: 0
success: True
njev: 4
nfev: 28
fun: -1.1623048291871221
x: array([ -3.60840218e-16, 2.24626781e-16, 1.63619563e-01, -2.27085639e-16, 8.36380437e-01])
message: 'Optimization terminated successfully.'
jac: array([ 1.81575805e-01, 5.40387481e-01, 8.18073750e-05, 1.03137662e+00, -1.60038471e-05, 0.00000000e+00])
nit: 4
得到的最优组合权重向量为：
In [10]:
opts['x'].round(3)
Out[10]:
array([-0. , 0. , 0.164, -0. , 0.836])
sharpe最大的组合3个统计数据分别为：
In [11]:
#预期收益率、预期波动率、最优夏普指数
statistics(opts['x']).round(3)
Out[11]:
array([ 0.508, 0.437, 1.162])
7.投资组合优化2——方差最小
接下来，我们通过方差最小来选出最优投资组合。
In [12]:
#但是我们定义一个函数对 方差进行最小化
def min_variance(weights):
return statistics(weights)[1]
optv = sco.minimize(min_variance, noa*[1./noa,],method = 'SLSQP', bounds = bnds, constraints = cons)
optv
Out[12]:
status: 0
success: True
njev: 7
nfev: 50
fun: 0.38542969450547221
x: array([ 1.14787640e-01, 3.28089742e-17, 2.09584008e-01, 3.53487044e-01, 3.22141307e-01])
message: 'Optimization terminated successfully.'
jac: array([ 0.3851725 , 0.43591119, 0.3861807 , 0.3849672 , 0.38553924, 0. ])
nit: 7
方差最小的最优组合权重向量及组合的统计数据分别为：
In [13]:
optv['x'].round(3)
Out[13]:
array([ 0.115, 0. , 0.21 , 0.353, 0.322])
In [14]:
#得到的预期收益率、波动率和夏普指数
statistics(optv['x']).round(3)
Out[14]:
array([ 0.226, 0.385, 0.587])
8.组合的有效前沿
有效前沿有既定的目标收益率下方差最小的投资组合构成。
在最优化时采用两个约束，1.给定目标收益率，2.投资组合权重和为1。
In [15]:
def min_variance(weights):
return statistics(weights)[1]
#在不同目标收益率水平（target_returns）循环时，最小化的一个约束条件会变化。
target_returns = np.linspace(0.0,0.5,50)
target_variance = []
for tar in target_returns:
cons = ({'type':'eq','fun':lambda x:statistics(x)[0]-tar},{'type':'eq','fun':lambda x:np.sum(x)-1})
res = sco.minimize(min_variance, noa*[1./noa,],method = 'SLSQP', bounds = bnds, constraints = cons)
target_variance.append(res['fun'])
target_variance = np.array(target_variance)
下面是最优化结果的展示。
叉号：构成的曲线是有效前沿（目标收益率下最优的投资组合）
红星：sharpe最大的投资组合
黄星：方差最小的投资组合
In [16]:
plt.figure(figsize = (8,4))
#圆圈：蒙特卡洛随机产生的组合分布
plt.scatter(port_variance, port_returns, c = port_returns/port_variance,marker = 'o')
#叉号：有效前沿
plt.scatter(target_variance,target_returns, c = target_returns/target_variance, marker = 'x')
#红星：标记最高sharpe组合
plt.plot(statistics(opts['x'])[1], statistics(opts['x'])[0], 'r*', markersize = 15.0)
#黄星：标记最小方差组合
plt.plot(statistics(optv['x'])[1], statistics(optv['x'])[0], 'y*', markersize = 15.0)
plt.grid(True)
plt.xlabel('expected volatility')
plt.ylabel('expected return')
plt.colorbar(label = 'Sharpe ratio')
