线性规划(Linear programming,简称LP),是运筹学中研究较早、发展较快、应用广泛、方法较成熟的一个重要分支,是辅助人们进行科学管理的一种数学方法,是研究线性约束条件下线性目标函数的极值问题的数学理论和方法。
线性规划是运筹学的一个重要分支,广泛应用于军事作战、经济分析、经营管理和工程技术等方面。为合理地利用有限的人力、物力、财力等资源作出的最优决策,提供科学的依据。
线性规划问题的解的概念
在人们的生产实践中,经常会遇到如何利用现有资源来安排生产,以取得最大经济效益的问题。此类问题构成了运筹学的一个重要分支—数学规划,而线性规划(Linear Programming 简记LP)则是数学规划的一个重要分支。自从1947 年G. B. Dantzig 提出求解线性规划的单纯形方法以来,线性规划在理论上趋向成熟,在实用中日益广泛与深入。特别是在计算机能处理成千上万个约束条件和决策变量的线性规划问题之后,线性规划的适用领域更为广泛了,已成为现代管理中经常采用的基本方法之一。
由于目标函数及约束条件均为线性函数,故被称为线性规划问题。线性规划问题是在一组线性约束条件的限制下,求一线性目标函数最大或最小的问题。
在数学上,线性规划问题的标准型为:
\max z=\sum^{n}_{j=1}c_jx_j\\ s.t.\left\{ \begin{aligned} &\sum^{n}_{j=1}a_{ij}x_j=b\ \ \ \ i=1,2,...,m\\ &x_j\geq0\ \ \ \ j=1,2,...,n\\ \end{aligned}\right.
其中上式被称为目标函数,下式被称为约束条件;
满足约束条件的解成为线性规划问题的可行解,而使目标函数达到最大值的可行解叫做最优解,所有可行解构成的集合称作问题的可行域,记为R
而在python或matlab程序中,规定的线性规划标准形式为:
\min_x\ c^T x \\ s.t.\left\{ \begin{aligned} &Ax \leq b \\ &Aeq·x = beq \\ &lb \leq x \leq ub \end{aligned}\right.
其中c和x为n维列向量,A、Aeq为适当维数的矩阵,b、 beq为适当维数的列向量。
Python实现线性规划的几种方式
Scipy库(简单,标准,常规)
Scipy 是一个用于数学、科学、工程领域的常用软件包,可以处理最优化、线性代数、积分、插值、拟合、特殊函数、快速傅里叶变换、信号处理、图像处理、常微分方程求解器等。
例题:求解下列线性规划问题
\max\ z=2x_1+3x_2-5x_3\\ s.t.\left\{ \begin{aligned} &x_1+x_2+x_3=7\\ &2x_1-5x_2+x_3\geq10\\ &x_1+3x_2+x_3\leq12\\ &x_1,x_2,x_3\geq0 \end{aligned}\right.
可以用Scipy库中的scipy.optimize(优化算法)模块进行求解线性规划问题,我们用到的函数为optimize.linprog(),其完整形式如下所示:
scipy.optimize.linprog(c,A_ub=None,b_ub=None,A_eq=None,b_eq=None,bounds=None,method='interior-point',callback=None,options=None,x0=None)
参数解释参见下表:
参数 | 解释 | 数据类型 | 备注 |
c | 线性目标函数的系数 | 一维数组 | 必选参数 |
A_ub | 不等式约束矩阵 | 二维数组 | |
b_ub | 不等式约束向量 | 一维数组 | |
A_eq | 等式约束矩阵 | 二维数组 | |
b_eq | 等式约束向量 | 一维数组 | |
bounds | x的范围,以(lb,ub)的形式表示 | 元组 | 1.使用None表示没有界限 2.默认情况下,界限为(0,None),即所有决策变量均为非负数 |
method | 算法,{‘interior-point’, ‘revised simplex’, ‘simplex’}三种算法可选 | 字符串 | 默认使用’interior-point’ |
callback | 调用回调函数 | ? | 不理解,建议别管 |
options | 求解器选项字典 | 字典 | 一般情况下用不到吧 |
x0 | 猜测决策变量的值,将通过优化算法进行优化。当前仅由’ revised simplex’ 方法使用此参数,并且仅当 x0 表示基本可行的解决方案时才可以使用此参数。 | 一维数组 | 一般情况下应该也用不到 |
一般来说函数会被重载为以下三种形式:
optimize.linprog(c,A_ub,b_ub)
optimize.linprog(c,A_ub,b_ub,Aeq,beq)
optimize.linprog(c,A_ub,b_ub,Aeq,beq,(LB,UB))
所以对于题给的限制条件,解法如下:
from scipy import optimize
import numpy as np
c = np.array([2,3,-5])
A_ub = np.array([[-2,5,-1],[1,3,1]])
b_ub = np.array([-10,12])
Aeq = np.array([[1,1,1]])
Beq = np.array([7])
#optimize.linprog(c,A,b,Aeq,beq,LB,UB,X0,OPTIOS)
res = optimize.linprog(-c,A_ub,b_ub,Aeq,Beq)
print(res)
返回的输出结果如下所示:
message: Optimization terminated successfully. (HiGHS Status 7: Optimal) success: True status: 0 fun: -14.571428571428571 x: [ 6.429e+00 5.714e-01 0.000e+00] nit: 3 lower: residual: [ 6.429e+00 5.714e-01 0.000e+00] marginals: [ 0.000e+00 0.000e+00 7.143e+00] upper: residual: [ inf inf inf] marginals: [ 0.000e+00 0.000e+00 0.000e+00] eqlin: residual: [ 0.000e+00] marginals: [-2.286e+00] ineqlin: residual: [ 0.000e+00 3.857e+00] marginals: [-1.429e-01 -0.000e+00] mip_node_count: 0 mip_dual_bound: 0.0 mip_gap: 0.0
对这些数据的解释如下:
输出内容 | 解释 | 数据类型 |
x | 在满足约束的情况下将目标函数最小化的决策变量的值 | 一维数组 |
fun | 目标函数的最佳值 | 浮点数 |
slack | 不等式约束的松弛值(名义上为正值)b_ub-A_ub*x | 一维数组 |
con | 等式约束的残差(名义上为零)b_eq-A_eq*x | 一维数组 |
success | 当算法成功找到最佳解决方案时为 True | 布尔值 |
status | 表示算法退出状态的整数0 : 优化成功终止1 : 达到了迭代限制2 : 问题似乎不可行3 : 问题似乎是不收敛4 : 遇到数值困难 | 整数 |
nit | 在所有阶段中执行的迭代总数 | 整数 |
message | 算法退出状态的字符串描述符 | 字符串 |
Pulp库(0-1规划问题,整数规划问题等)
Python的scipy库中提供了解简单线性或非线性规划问题,但是不能求解如背包问题的0-1规划问题,或整数规划问题,混合整数规划问题,pulp库可以求解以上类型的问题,并且有更多的通用性,编写程序更自由。
pulp在使用上 首先是定义一个问题,第一个参数为问题的名称,不过并非必要的参数,而通过sense参数可决定优化的是最大值还是最小值
prob = pulp.LpProblem('problem name', sense=pulp.LpMinimize)
然后是定义变量:
x0 = pulp.LpVariable('x0', lowBound=0, upBound=None, cat=pulp.LpInteger)
x1 = pulp.LpVariable('x1', lowBound=0, upBound=None, cat=pulp.LpInteger)
x2 = pulp.LpVariable('x2', lowBound=0, upBound=None, cat=pulp.LpInteger)
这里定义了三个变量,第一个参数为变量名,lowBound 、upBound为上下限,cat为变量类型,默认为连续变量,还可以设为离散变量或二值变量,具体怎么设置?
如上述代码所示,pulp.LpInteger为离散变量,pulp.LpBinary为二值变量,同时也可以传入’Integer’字符串的方式指明变量类型。从上面几个问题的代码可以看出,我几乎没有单个定义变量,而是批量定义的。
然后是添加目标函数:
prob += 2*x0-5*x1+4*x2
只要让问题(prob)直接加变量的表达式即可添加目标函数。
再之后是添加约束:
prob += (x0+x1-6*x2 <= 120)
只要让问题(prob)直接加变量的判断式即可添加约束
调用solve方法解出答案,如果省去这一步,后面的变量和结果值都会显示None。
prob.solve()
print(pulp.value(prob.objective))
print(pulp.value(x0))
打印优化结果,并显示当优化达到结果时x0的取值。
例题:求解下列线性规划问题
\max\ z=2x_1+3x_2-5x_3\\ s.t.\left\{ \begin{aligned} &x_1+x_2+x_3=7\\ &2x_1-5x_2+x_3\geq10\\ &x_1+3x_2+x_3\leq12\\ &x_1,x_2,x_3\geq0 \end{aligned}\right.
问题和scipy库一样,但是用pulp库来求解,代码如下:
import pulp
#目标函数的系数
z=[2,3,-5]
#约束
a=[[2,-5,1],[-1,-3,-1]]
b=[10,-12]
A_eq=[1,1,1]
b_eq=7
#确定最大化最小化问题,最小化只要把Max改成Min即可
m = pulp.LpProblem(sense=pulp.LpMaximize)
#定义三个变量放到列表中
x = [pulp.LpVariable(f'x{i}',lowBound=0) for i in [1,2,3]]
#定义目标函数,lpDot可以将两个列表的对应位相乘再加和
#相当于z[0]*x[0]+z[1]*x[0]+z[2]*x[2]
m += pulp.lpDot(z,x)
#设置约束条件
for i in range(len(a)):
m+=(pulp.lpDot(a[i],x)>=b[i])
m+=(pulp.lpDot(A_eq,x)==b_eq)
#求解
m.solve()
#输出结果
print(f'优化结果:{pulp.value(m.objective)}')
print(f'参数取值:{[pulp.value(var) for var in x]}')
参考的输出结果如下所示:
优化结果:14.57142851 参数取值:[6.4285714, 0.57142857, 0.0]
Cvxpy库(凸优化问题可用)
CVXPY是一种可以内置于Python中的模型编程语言,解决凸优化问题(整数规划、01规划和混合规划)。它可以自动转化问题为标准形式,调用解法器,解包结果集。cvxpy包相对前面2种算是最专业了,功能也更强大。
例题:求解下列线性规划问题:
\min z=13x_1+9x_2+10x_3+11x_4+12x_5+8x_6\\ s.t.\left\{ \begin{aligned} &x_1+x_4=400\\ &x_2+x_5=600\\ &x_3+x_6=500\\ &0.4x_1+1.1x_2+x_3 \leq 800\\ &0.5x_4+1.2x_5+1.3x_6 \leq 900\\ &x_i \geq 0,i=1,2,...,6 \end{aligned}\right.
cvxpy库解决线性规划问题的思路略有不同,具体的实现方式可以参照以下代码:
import cvxpy as cp
import numpy as np
# 定义自变量
n = 6 # 有6个自变量:x1,x2,x3,x4,x5,x6
x = cp.Variable(n)
# 定义相等约束条件
A1 = np.array([[1, 0, 0, 1, 0, 0],
[0, 1, 0, 0, 1, 0],
[0, 0, 1, 0, 0, 1],
])
b1 = np.array([400, 600, 500])
#定义小于等于约束条件
A2 = np.array([[0.4, 1.1, 1, 0, 0, 0],
[0, 0, 0, 0.5, 1.2, 1.3],
])
b2 = np.array([800, 900])
#定义大于等于约束条件
A3 = np.eye(n)
b3 = np.zeros(n)
#定义目标函数
c = np.array([13, 9, 10, 11, 12, 8])
# 定义问题,把三个约束条件添加上
prob = cp.Problem(cp.Minimize(c.T @ x),
[A1 @ x == b1, A2 @ x <= b2, A3 @ x >= b3])
# 解决问题
prob.solve()
# 输出结果
print("\n目标函数的最小值", prob.value)
print("所求的x矩阵")
print(x.value)
# 对x向量各元素取整数后再输出
for item in x.value:
print(round(item))
参考的输出结果如下所示:
目标函数的最小值 13800.00000202733 所求的x矩阵 [1.85735461e-07 6.00000000e+02 1.50597739e-07 4.00000000e+02 4.51551936e-07 5.00000000e+02] 0 600 0 400 0 500
一些可以转化为线性规划的问题
目标函数含绝对值的情况
例1.求解下列线性规划问题
\min\ |x_1|+|x_2|+...+|x_n|\\ s.t.\ Ax\leq b\\ 其中x=[x_1\ ...\ x_n]^T,A和b为相应维数的矩阵和向量
可以采取这样的思路进行转化:
\begin{aligned} &取u_i=\frac{x_i+|x_i|}{2},v_i=\frac{|x_i|-x_i}{2},记u=[u_1\ ...\ u_n]^T,v=[v_1\ ...\ v_n]^T\\ &可以将问题转换为如下线性规划模型: \end{aligned}\\ \min \sum^{n}_{i=1}(u_i+v_i)\\ s.t.\left\{ \begin{aligned} &A(u-v) \leq b\\ &u,v \geq 0 \end{aligned}\right.
如果用代码来实现的话,大致是这样的(数据可以根据具体情况修改):
from scipy import optimize
import numpy as np
x=np.array([1,2,3,4,5])
u=(abs(x)+x)/2
v=(abs(x)-x)/2
c=u-v
A_ub = np.array([[-2,5,-1,5,1],[1,3,1,4,3]])
b_ub = np.array([-10,12])
res = optimize.linprog(-c,A_ub,b_ub)
print(res)
运输问题
某商品有m个产地、n个销售地,各产地的产量分别为a1,a2,…,am各销售地的需求量分别为b1,…,bn。若该商品由i产地运输到j销售地的单位运输价格为c_{ij},问应该如何调运才能使总运费最省?
引入变量x_{ij},其取值为由i产地运往j销售地的该商品数量,数学模型为
\min \sum^{m}_{i=1} \sum^{n}_{j=1}c_{ij}x_{ij}\\ s.t.\left\{ \begin{aligned} &\sum^{n}_{j=1}x_{ij}=a_i,i=1,2,...,m\\ &\sum^{m}_{i=1}x_{ij}=b_j,j=1,2,...,n\\ \end{aligned}\right.
例2:一个农民承包了6块耕地共300亩,准备播种小麦、玉米、水果和蔬菜四种农产品,各种农产品的计划播种面积、每块土地种植不同农产品的不同收益如下表:
地块1 | 地块2 | 地块3 | 地块4 | 地块5 | 地块6 | 计划播种面积(亩) | |
小麦 | 500 | 550 | 630 | 1000 | 800 | 700 | 76 |
玉米 | 800 | 700 | 600 | 950 | 900 | 930 | 88 |
水果 | 1000 | 960 | 840 | 650 | 600 | 700 | 96 |
蔬菜 | 1200 | 1040 | 980 | 860 | 880 | 780 | 40 |
地块面积(亩) | 42 | 56 | 44 | 39 | 60 | 59 |
问如何安排种植计划,可以得到最大的总收益
用pulp库解决这一问题的代码如下所示:
import pulp
import numpy as np
from pprint import pprint
#运输问题解决函数
def transportation_problem(costs, x_max, y_max):
row = len(costs)
col = len(costs[0])
prob = pulp.LpProblem('Transportation Problem', sense=pulp.LpMaximize)
var = [[pulp.LpVariable(f'x{i}{j}', lowBound=0, cat=pulp.LpInteger) for j in range(col)] for i in range(row)]
flatten = lambda x: [y for l in x for y in flatten(l)] if type(x) is list else [x]
prob += pulp.lpDot(flatten(var), costs.flatten())
for i in range(row):
prob += (pulp.lpSum(var[i]) <= x_max[i])
for j in range(col):
prob += (pulp.lpSum([var[i][j] for i in range(row)]) <= y_max[j])
prob.solve()
return {'objective':pulp.value(prob.objective), 'var': [[pulp.value(var[i][j]) for j in range(col)] for i in range(row)]}
#传入数据
if __name__ == '__main__':
costs = np.array([[500, 550, 630, 1000, 800, 700],
[800, 700, 600, 950, 900, 930],
[1000, 960, 840, 650, 600, 700],
[1200, 1040, 980, 860, 880, 780]])
max_plant = [76, 88, 96, 40]
max_cultivation = [42, 56, 44, 39, 60, 59]
res = transportation_problem(costs, max_plant, max_cultivation)
print(f'最大值为{res["objective"]}')
print('各变量的取值为:')
pprint(res['var'])
指派问题(1-0规划问题)
拟分配n人去干n项工作,没人干且仅干一项工作,若分配第i人去干第j项工作,需要花费c_ij单位时间,问应该如何分配工作才能使公认花费的总时间最少?
引入变量x_ij,若分配i干j工作,则取x_ij=1,否则取0,则该指派问题可以简化成模型为
\min \sum^{n}_{i=1} \sum^{n}_{j=1}c_{ij}x_{ij}\\ s.t.\left\{ \begin{aligned} &\sum^{n}_{j=1}x_{ij}=1\\ &\sum^{n}_{i=1}x_{ij}=1\\ &0 \leq x_{ij} \leq 1\\ \end{aligned}\right.
例3.假设指派问题的系数矩阵如下所示,解决该指派问题
C=\left [ \begin{aligned} &12\quad&&7\quad&&9\quad&&7\quad&&9\\ &8\quad&&9\quad&&6\quad&&6\quad&&6\\ &12\quad&&17\quad&&12\quad&&14\quad&&12\\ &12\quad&&14\quad&&6\quad&&16\quad&&10\\ &12\quad&&10\quad&&7\quad&&0\quad&&6\\ \end{aligned}\right]
Scipy库的解决方法
代码如下:
#引入包,linear_sum_assignment是专门用于解决指派问题的模块
import numpy as np
from scipy.optimize import linear_sum_assignment
#导入系数矩阵
efficiency_matrix = np.array([
[12,7,9,7,9],
[8,9,6,6,6],
[7,17,12,14,12],
[15,14,6,6,10],
[4,10,7,10,6]])
#定义了开销矩阵(指派问题的系数矩阵)efficiency_matrix,传入linear_sum_assignment,
#结果返回的是最优指派的行和列,例如第一行选择第二列,意为:将第一个人派往第二个工作。
#而根据numpy.array的性质,传入行和列就会返回行列所对应的值,即为输出的第三列
row_index, col_index=linear_sum_assignment(efficiency_matrix)
print(row_index+1)
print(col_index+1)
print(efficiency_matrix[row_index,col_index])
#对其求和,即可得到指派问题的最小时间
print(efficiency_matrix[row_index, col_index].sum())
参考输出结果如下:
[1 2 3 4 5] [2 3 1 4 5] [7 6 7 6 6] 32
Pulp库的解决方法
代码如下:
import pulp
import numpy as np
#先定义通用解决方法,其中的flatten是递归展开列表用的
def assignment_problem(efficiency_matrix):
row = len(efficiency_matrix)
col = len(efficiency_matrix[0])
flatten = lambda x: [y for l in x for y in flatten(l)] if type(x) is list else [x]
m = pulp.LpProblem('assignment', sense=pulp.LpMinimize)
var_x = [[pulp.LpVariable(f'x{i}{j}', cat=pulp.LpBinary) for j in range(col)] for i in range(row)]
m += pulp.lpDot(efficiency_matrix.flatten(), flatten(var_x))
for i in range(row):
m += (pulp.lpDot(var_x[i], [1]*col) == 1)
for j in range(col):
m += (pulp.lpDot([var_x[i][j] for i in range(row)], [1]*row) == 1)
m.solve()
return {'objective':pulp.value(m.objective), 'var': [[pulp.value(var_x[i][j]) for j in range(col)] for i in range(row)]}
#然后定义矩阵,输入求解
efficiency_matrix = np.array([
[12, 7, 9, 7, 9],
[8, 9, 6, 6, 6],
[7, 17, 12, 14, 9],
[15, 14, 6, 6, 10],
[4, 10, 7, 10, 9]
])
res = assignment_problem(efficiency_matrix)
print(f'最小值{res["objective"]}')
print(res['var'])
参考输出结果如下:
最小值32.0 [[0.0, 1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 1.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0, 0.0]]
线性规划综合实战:投资的收益和风险模型
问题提出
市场上有n种资产s_i(i=1,2,…,n)可以选择,现用数额为M的相当大的资金作一个时期的投资。这n种资产在这一时期内购买s_i的平均收益率为r_i,风险损失率为q_i;,投资越分散,总的风险越少,总体风险可用投资的s;中最大的一个风险来度量。
购买s_i时要付交易费,费率为p_i;,当购买额不超过给定值u_i时,交易费按购买u_i计算。另外,假定同期银行存款利率是r_0%,既无交易费又无风险(r_0=5%)。
已知n=4时相关数据如表1.1。
s_i | r_i(%) | q_i(%) | p_i(%) | u_i(元) |
s_1 | 28 | 2.5 | 1 | 103 |
s_2 | 21 | 1.5 | 2 | 198 |
s_3 | 23 | 5.5 | 4.5 | 52 |
s_4 | 25 | 2.6 | 6.5 | 40 |
试给该公司设计一种投资组合方案,即用给定资金M,有选择地购买若干种资产或存银行生息,使净收益尽可能大,使总体风险尽可能小。
符号的基本规定与假设
符号规定
s_i表示第i种投资项目,如股票,债券等,i= 0,1,… ,n,其中s_0指存入银行;
r_i, p_i,q_i,分别表示s_i的平均收益率,交易费率,风险损失率,i= 0,…,n,
其中p_0=0,q_0=0;
u_i表示s_i的交易定额,i= 1,…,n;
x_i表示投资项目s_i的资金,i= 0,1,…,n ;
a表示投资风险度;
Q表示总体收益;
基本假设
(1)投资数额M相当大,为了便于计算,假设M=1;
(2)投资越分散,总的风险越小;
(3)总体风险用投资项目s_i中最大的一个风险来度量;
(4) n+1种资产s_i之间是相互独立的;
(5)在投资的这一时期内,r_i,p_i,q_i,为定值,不受意外因素影响;
(6)净收益和总体风险只受r_i,p_i,q_i影响,不受其它因素干扰。
模型的分析与建立
1.总体风险用所投资的s_i中最大的一个风险来衡量,即
\max\{q_ix_i\ |\ i=1,2,...,n\}\\
2.购买s_i(i=1,…,n)所付交易费是一个分段函数,即
交易费=\left\{\begin{aligned} &p_ix_i,\ \ \ \ x_i>u_i\\ &p_iu_i,\ \ \ \ x_i \leq u_i \end{aligned}\right.\\
3.要使净收益尽可能大,总体风险尽可能小,这是一个多目标规划模型
目标函数为
\left\{\begin{aligned} &\max \sum^{n}_{i=0}(r_i-p_i)x_i\\ &\min \max \{q_ix_i\}\\ \end{aligned}\right.\\
约束条件为
\left\{\begin{aligned} &\sum^{n}_{i=0}(1+p_i)x_i=M\\ &x_i \geq 0,\ \ i=0,1,...,n\\ \end{aligned}\right.\\
下文介绍了几种方法将该多目标规划转化为单一目标的线性规划,但是只引用对该模型的介绍以及附上我自己对模型的求解代码,不对求解结果作具体分析
模型一 固定风险水平,优化收益
在实际投资中,投资者承受风险的程度不一样,若给定风险一个界限a,使得最大的一个风险(q_i*x_i)/M<=a,可找到相应的投资方案。这样把多目标规划变成一个目标的线性规划。
该模型可以表示为以下这种形式:
\max \sum^{n}_{i=0}(r_i-p_i)x_i\\ s.t.\left\{\begin{aligned} & \frac{q_ix_i}{M} \leq a\\ &\sum^{n}_{i=0}(1+p_i)x_i=M,\ \ x_i\geq0,\ \ i=0,1,...,n \end{aligned}\right.\\
用python的实现方法如下所示:
from scipy import optimize
import numpy as np
from matplotlib import pyplot as plt
#题目给到的数据
r=0.01*np.array([5,28,21,23,25]) #平均收益率
q=0.01*np.array([0,2.5,1.5,5.5,2.6]) #交易费率
p=0.01*np.array([0,1,2,4.5,6.5]) #风险损失率
u=np.array([103,198,52,40]) #s_i的交易定额
M=1 #根据基本假设为了方便计算使M=1
a=0 #投资风险度是由用户给出的变量,暂时设初值为0
Qs=np.array([]) #总收益的最优解集合,随a的变化而变化
#对(0,0.05)之间的a进行遍历操作
while a<0.05:
c=r-p
A_ub=q*np.eye(len(c))
b_ub=a*np.ones(len(c))
Aeq=np.array([1+p])
beq=np.array([1])
Q_i=optimize.linprog(-c,A_ub,b_ub,Aeq,beq)
Qs=np.append(Qs,-Q_i.fun)
a+=0.001
#绘图,得出结果
plt.title("Risk & Reward")
plt.xlabel("a")
plt.ylabel("Q")
plt.plot(np.arange(0,0.05,0.001),Qs,"ob")
plt.show()
模型二 固定盈利水平,极小化风险
若投资者希望总盈利至少达到水平k以上,在风险最小的情况下寻求相应的投资组合。
该模型可以表示为以下这种形式:
\min\{\max\{q_ix_i\}\}\\ s.t.\left\{\begin{aligned} & \sum^{n}_{i=0}(r_i-p_i)x_i \geq k\\ &\sum^{n}_{i=0}(1+p_i)x_i=M,\ \ x_i\geq0,\ \ i=0,1,...,n \end{aligned}\right.\\
由于目标函数使用了max函数,使用scipy和pulp库的解决方法我想不太到,所以考虑使用cvxpy库来解决该问题,代码如下:
from matplotlib import pyplot as plt
import numpy as np
import cvxpy as cp
#题目给到的数据
r=0.01*np.array([5,28,21,23,25]) #平均收益率
q=0.01*np.array([0,2.5,1.5,5.5,2.6]) #交易费率
p=0.01*np.array([0,1,2,4.5,6.5]) #风险损失率
u=np.array([103,198,52,40]) #s_i的交易定额
M=1 #根据基本假设为了方便计算使M=1
k=0.05 #期望总收益是由用户给出的变量,暂时设初值为0.05
kk=[] #k的集合,自变量轴
Rs=[] #总体风险的最优解集合,随k的变化而变化
x=cp.Variable(6,pos=True) #决策向量6个,且都为正数
obj=cp.Minimize(x[5]) #优化得目标函数(构造目标函数)
#对(0.05,0.3)之间的k进行遍历操作
while k<0.3:
con=[cp.multiply(q[1:5],x[1:5])-x[5]<=0,(r-p)@x[:-1]>=k,(1+p)@x[:-1]==1]
prob=cp.Problem(obj,con)
prob.solve()
kk.append(k)
Rs.append(prob.value)
k=k+0.005
#绘图,得出结果
plt.figure(figsize=(20,8),dpi=80)
plt.title("Risk & Reward")
plt.xlabel("k")
plt.ylabel("R")
plt.plot(kk,Rs,"ob")
plt.show()
大概可以看出来这段代码实现max函数的方式是再引入一个决策变量x[5]作为q*x[0:4]的最大值,最大值体现在额外的附加条件cp.multiply(q[1:5],x[1:5])-x[5]<=0,不过说实话我觉得这里应该是x[0:4]才对,但是最后结果证明x[1:5]才是正解,所以就别去试图理解它了吧)
模型三 投资偏好,同时考虑风险与收益
投资者在权衡资产风险和预期收益两方面时,希望选择一个令自己满意的投资组合。因此对风险、收益分别赋予权重s(0<s≤1)和(1- s),s称为投资偏好系数。
该模型可以表示为以下这种形式:
\min\ s\{\max\{q_ix_i\}\}-(1-s)\sum^{n}_{i=0}(r_i-p_i)x_i\\ s.t.\ \sum^{n}_{i=0}(1+p_i)x_i=M,\ \ x_i\geq0,\ \ i=0,1,...,n
代码我暂时性开摆了,等有空再来填坑))