用于优化建模的Python代数工具箱
pyatom的Python项目详细描述
PyAtom公司
PyAtom是一个开源的Python代数工具箱，用于优化建模。它与NumPy包在N维数组操作方面是一致的，并且依赖于商业求解器Gurobi来解决格式化模型。到目前为止，该工具箱能够对格式化为二阶锥程序的鲁棒优化问题进行建模。在
示例
线性规划
我们使用一个非常简单的线性程序来演示运行PyAtom来建模优化问题的一般过程。在
importpyatom.lpaslp# Import the LP modeling module from PyAtomimportpyatom.grb_solverasgrb# Import the solver interface for Gurobimodel=lp.Model()# Create an LP modelx=model.dvar()# Define a decision variable xy=model.dvar()# Define a decision variable ymodel.max(3*x+4*y)# Maximize the objective functionmodel.st(2.5*x+y<=20)# Define constraintsmodel.st(3*x+3*y<=30)model.st(x+2*y<=16)model.st(x<=3)model.st(abs(y)<=2)model.solve(grb)# Solve the model by Gurobi^{pr2}$
通过模型对象或变量对象的get（）方法，可以得到最优目标值和最优解。在
print(model.get())print(x.get())print(y.get())
17.0
[3.]
[2.]
model.do_math()
Linear program object:
=============================================
Number of variables: 3
Continuous/binaries/integers: 3/0/0
---------------------------------------------
Number of linear constraints: 6
Inequalities/equalities: 6/0
Number of coefficients: 11
请注意，决策变量的最优解是以数组类型的对象给出的。为了便于调试，PyAtom包可以通过运行命令生成一个约束信息表数学模型（）.showlc（）。在
model.do_math().showlc()
| x1 | x2 | x3 | sense | constants | |
|---|---|---|---|---|---|
| LC1 | 0.0 | 2.5 | 1.0 | <= | 20.0 |
| LC2 | 0.0 | 3.0 | 3.0 | <= | 30.0 |
| LC3 | 0.0 | 1.0 | 2.0 | <= | 16.0 |
| LC4 | -1.0 | 3.0 | -4.0 | <= | 0.0 |
| LC5 | 0.0 | 0.0 | 1.0 | <= | 2.0 |
| LC6 | 0.0 | 0.0 | -1.0 | <= | 2.0 |
均值-方差投资组合优化
在这个例子中，我们将均值-方差投资组合优化问题描述为一个二阶锥规划。这个模型的细节可以从论文Price of Robustness中找到。在
importpyatomasat# Import the PyAtom packageimportpyatom.socpascp# Import the SOCP modeling module from PyAtomimportpyatom.grb_solverasgrb# Import the solver interface for Gurobiimportnumpyasnpn=150# Number of stocksi=np.arange(1,n+1)p=1.15+0.05/n*i# Mean values of stock returnssig=0.05/450*(2*i*n*(n+1))**0.5# Standard deviation of stock returnsphi=5# Constant phimodel=cp.Model()# Create an SOCP modelx=model.dvar(n)# Define decision variables as an array with length nmodel.max(p@x-phi*at.sumsqr(sig*x))# Define the objective functionmodel.st(sum(x)==1)# Define the constraintsmodel.st(x>=0)model.solve(grb)# Solve the model by Gurobi
Being solved by Gurobi...
Solution status: 2
Running time: 0.0024s
关于最优库存分配的解决方案如下所示。在
importmatplotlib.pyplotaspltplt.plot(i,x.get())plt.xlabel('stocks')plt.ylabel('x')plt.show()
稳健投资组合优化
文中介绍的稳健投资组合优化模型Price of Robustness也可以用PyAtom包来表示。在
importpyatom.roasro# Import the robust optimization module from PyAtomn=150# Number of stocksi=np.arange(1,n+1)p=1.15+0.05/n*i# Mean values of stock returnssig=0.05/450*(2*i*n*(n+1))**0.5# Maximum deviation of stock returnsGamma=3# Budget of uncertaintymodel=ro.Model()# Create a robust optimization modelx=model.dvar(n)# Define decision variables x as an arrayz=model.rvar(n)# Define random variables z as an arraymodel.maxmin((p+sig*z)@x,# Define the max-min objective functionabs(z)<=1,at.norm(z,1)<=Gamma)# Uncertainty set of the modelmodel.st(sum(x)<=1)# Define constraintsmodel.st(x>=0)model.solve(grb)# Solve the model by Gurobi
Being solved by Gurobi...
Solution status: 2
Running time: 0.0029s
股票分配的解决方案如下所示。在
plt.plot(i,x.get())plt.xlabel('stocks')plt.ylabel('x')plt.show()
背包问题：鲁棒优化v.s.鲁棒优化
在这个例子中，我们使用PyAtom包来实现本文The Dao of Robustness中描述的健壮性和健壮性优化模型。在
importpyatom.roasroimportpyatom.grb_solverasgrbimportpyatomasatimportnumpyasnpimportnumpy.randomasrdimportmatplotlib.pyplotaspltN=50b=2000c=2*rd.randint(low=5,high=10,size=N)# Profit coefficients ofw_hat=2*rd.randint(low=10,high=41,size=N)# Baseline values of the weightsdelta=0.2*w_hat# Maximum deviations
鲁棒优化方法的函数如下所示。在
defrobust(r):""" The function robust implements the robust optmization model, given the budget of uncertainty r """model=ro.Model('robust')x=model.dvar(N,vtype='B')# Define decision variablesz=model.rvar(N)# Define random variablesmodel.max(c@x)model.st(((w_hat+z*delta)@x<=b).forall(abs(z)<=1,at.norm(z,1)<=r))model.solve(grb,display=False)returnmodel.get(),x.get()# Return objective value and the optimal solution
鲁棒优化模型的函数。在
defrobustness(target):model=ro.Model('robustness')x=model.dvar(N,vtype='B')k=model.dvar()z=model.rvar(N)u=model.rvar(N)model.min(k)model.st(c@x>=target)model.st(((w_hat+z*delta)@x-b<=k*u.sum()).forall(abs(z)<=u,u<=1))model.st(k>=0)model.solve(grb,display=False)returnmodel.get(),x.get()
以下函数sim用于通过模拟计算违规概率。在
defsim(x_sol,zs):""" The function sim is for calculating the probability of violation via simulations. x_sol: solution of the Knapsack problem zs: random sample of the random variable z """ws=w_hat+zs*deltareturn(ws@x_sol>b).mean()
通过使用上述函数，我们运行鲁棒和鲁棒优化模型。在
step=0.1rs=np.arange(1,5+step,step)# All budgets of uncertaintynum_samp=20000zs=1-2*rd.rand(num_samp,N)# Random samples for z"""Robust optimization"""outputs_rb=[robust(r)forrinrs]# Robust optimization modelstgts=[output[0]foroutputinoutputs_rb]# Objective values as the targetspv_rb=[sim(output[1],zs)foroutputinoutputs_rb]# Prob. violation for robust optimization"""Robustness optimization"""outputs_rbn=[robustness(tgt)fortgtintgts]# Robustness optimization modelspv_rbn=[sim(output[1],zs)foroutputinoutputs_rbn]# Objective values as the targets
结果显示如下。在
plt.plot(rs,pv_rb,marker='o',markersize=5,c='b',label='Robust Optimization')plt.plot(rs,pv_rbn,c='r',label='Robustness Optimization')plt.legend()plt.xlabel('Parameter r in robust optimization')plt.ylabel('Prob. violation')plt.show()plt.scatter(tgts,pv_rb,c='b',alpha=0.3,label='Robust Optimization')plt.scatter(tgts,pv_rbn,c='r',alpha=0.3,label='Robustness Optimization')plt.legend()plt.xlabel(r'Target return $\tau$')plt.ylabel('Prob. violation')plt.show()
线性决策规则的鲁棒优化
frompyatomimportroimportpyatom.grb_solverasgrbimportnumpyasnpimportnumpy.randomasrdimportmatplotlib.pyplotaspltnp.set_printoptions(precision=2,linewidth=200)
批量问题的鲁棒优化
用PyAtom实现了本文中描述的批量问题的鲁棒模型，以演示如何为自适应决策指定线性决策规则。在
N=30情况下的参数定义如下。在
^{pr21}$然后由以下代码实现该模型。在
model=ro.Model()d=model.rvar(N)x=model.dvar(N)y=model.ldr((N,N))y.adapt(d)model.minmax((c*x).sum()+(tmat*y).sum(),(d>=0,d<=dmax,sum(d)<=Gamma))model.st(d<=y.sum(axis=0)-y.sum(axis=1)+x)model.st(y>=0)model.st(x>=0)model.st(x<=20)model.solve(grb)
Being solved by Gurobi...
Solution status: 2
Running time: 0.4217s
可以通过调用方法x.get（）来检索股票分配x的最优解，结果如下所示。在
plt.figure(figsize=(5,5))plt.scatter(xy[0],xy[1],c='w',edgecolors='k')plt.scatter(xy[0],xy[1],s=50*x.get(),c='k',alpha=0.4)plt.axis('equal')plt.xlim([-1,11])plt.ylim([-1,11])plt.grid()plt.show()
库存
本文用本文中的AARC库存模型Adjustable robust solutions of uncertain linear programs来演示多阶段问题线性决策规则的实现。在
frompyatomimportroimportnumpyasnpimportpyatomasatimportpyatom.grb_solverasgrbT,I=24,3P,Q=567,13600V0,Vmin,Vmax=1500,500,2000delta=0.2t=np.arange(T).reshape((1,T))d0=1000*(1+1/2*np.sin(np.pi*(t/12)))c=np.array([[1],[1.5],[2]])@(1+1/2*np.sin(np.pi*t/12))model=ro.Model()d=model.rvar(T)# Random demanduset=(d<=(1+delta)*d0,d>=(1-delta)*d0)# Uncertainty set of random demandp=model.ldr((I,T))# decision rule for adaptive decisionforiinrange(1,T):p[:,i].adapt(d[:i])# Define LDR dependencymodel.minmax((c*p).sum(),uset)# The min-max objective functionmodel.st(p>=0);# Constraintsmodel.st(p<=P);model.st(p.sum(axis=1)<=Q)v=V0foriinrange(T):v+=p[:,i].sum()-d[i]model.st(v<=Vmax)model.st(v>=Vmin)model.solve(grb)
Being solved by Gurobi...
Solution status: 2
Running time: 2.0997s
- 项目
标签：
欢迎加入QQ群-->： 979659372
