4. 最优化建模

• 理论上可以用最优性条件求“非线性规划”的最优解，但在实践中并不切实可行。
• 一般情况下，求解最优化问题所用的计算方法最常见的是“迭代下降算法”。

最优性条件：问题的最优解所满足的必要或者充分条件。它为最优化问题求解算法的设计、分析提供必不可少的理论基础。

无约束问题的极值条件

• 一阶必要条件：设目标函数在点处可微，若是局部极小点，则.
• 二阶必要条件：设目标函数在点处二次可微，若是局部极小点，则, 并且Hessian矩阵.
• 二阶充分条件：设目标函数在点处二次可微，若且, 则是局部极小点。
• 充要条件：设是定义在上的可微凸函数，则是整体极小点（全局最优解）的充要条件是.

约束问题的最优性条件

Kuhn-Tucker必要条件：设为约束问题的可行点， 记, 和在点可微， 在点连续，在点连续可微， 向量集线性无关。如果是局部最优解， 则存在数和使得

定义Lagrange函数

若为问题局部最优解，则存在乘子向量使得 .

此时，一阶必要条件可以表达为

求：

转换为标准形

import numpy as np

z = np.array([4, 3])
a = np.array([[2, 1.0], [1,1], [0,1]])
b = np.array([10,8, 7])
x1_bound = x2_bound =(0, None)

from scipy import optimize

res = optimize.linprog(-z, A_ub=a, b_ub=b,bounds=(x1_bound, x2_bound))
print(res)

结果

     fun: -26.0
 message: 'Optimization terminated successfully.'
     nit: 2
   slack: array([ 0.,  0.,  1.])
  status: 0
 success: True
       x: array([ 2.,  6.])

假设

1. 为便于计算，设投资总额为。
2. 在整个投资时期内，, , 为定值，不受意外因素的影响。
3. 净收入和整体风险只受, , 的影响
4. 假定总投资额很大，每个产品投资总额会远大于，这样，交易费简单地用来计算。
5. 总风险按投资产品中，风险最大的来算。
6. 银行存款计为产品

使收益尽可能大，风险尽可能小，这是个多目标规划模型。

目标函数：

约束条件：

模型一 固定风险水平。实际投资中，给出投资者的风险界限，也就是说，最多可以承受的资金损失，即。因此，模型为

模型二 固定盈利水平，极小化风险。实际投资中，投资者希望总盈利达到，找一个风险最小的组合。模型为

模型三 在风险和收益间，找一个平衡点。对风险、收益分别赋予权重和。模型为

解模型一

针对固定风险水平从到，每获取一个结果，

可以看到，在到，收益增加最快

0.0 0.05
0.001 0.0755276107226
0.002 0.101055221445
0.003 0.126582832168
0.004 0.15211044289
0.005 0.177638053613
0.006 0.201907639778
0.007 0.206607426376
0.008 0.21124338118
0.009000000000000001 0.215519617225
0.010000000000000002 0.219019607843
0.011000000000000003 0.222294117647
0.012000000000000004 0.225568627451
0.013000000000000005 0.228843137255
0.014000000000000005 0.232117647059

模型

设种玉米亩，种小麦亩，种燕麦亩，则

目标：

约束：

Mathematica

Maximize[{400 x1 + 200 x2 + 150 x3, x1 + x2 + x3 <= 625, 
  0.8 x1 + 0.2 x2 + 0.3 x3 <= 300, 3.0 x1 + 1.0 x2 + 1.5 x3 <= 1000, 
  x1 >= 0, x2 >= 0, x3 >= 0}, {x1, x2, x3}]

得到结果

{162500., {x1 -> 187.5, x2 -> 437.5, x3 -> 0.}}

模型

设120亩地中，种玉米亩，种小麦亩，种燕麦亩；25亩地中，种玉米亩，种小麦亩，种燕麦亩；则

目标：

约束：

Mathematica

Maximize[{120 (400 x1 + 200 x2 + 150 x3) + 
   25 (400 x4 + 200 x5 + 150 x6), x1 + x2 + x3 <= 5, 
  x4 + x5 + x6 <= 1, 
  120 (0.8 x1 + 0.2 x2 + 0.3 x3) + 25 (0.8 x4 + 0.2 x5 + 0.3 x6) <= 
   300, 120 (3.0 x1 + 1.0 x2 + 1.5 x3) + 
    25 (3.0 x4 + 1.0 x5 + 1.5 x6) <= 1000, x1 >= 0, x2 >= 0, x3 >= 0, 
  x4 >= 0, x5 >= 0, x6 >= 0, 
  Element[{x1, x2, x3, x4, x5, x6}, Integers]}, {x1, x2, x3, x4, x5, 
  x6}]

得到结果

{154000., {x1 -> 1, x2 -> 4, x3 -> 0, x4 -> 1, x5 -> 0, x6 -> 0}}

模型

设生产小型车辆，中型车车，大型车辆，则

目标:

约束:

方法一： 将问题分解为8个线性规划子模型

逐一求解这8个问题，取出最优解

方法二：变为非线性规划。 把约束条件变为

引入0-1变量： 设, , 只取或两个值，则条件变为

其中为充分大的数。

Mathematica

In[279]:= Maximize[{2 x1 + 3 x2 + 4 x3, 1.5 x1 + 3 x2 + 5 x3 <= 600, 
  280 x1 + 250 x2 + 400 x3 <= 60000, x1 <= 10^9*y1, x2 <= 10^9 y2, 
  x3 <= 10^9  y3, x1 >= 80 y1, x2 >= 80 y2, x3 >= 80 y3, x1 >= 0, 
  x2 >= 0, x3 >= 0, Element[{x1, x2, x3}, Integers], 
  y1 == 0 || y1 == 1, y2 == 0 || y2 == 1}, {x1, x2, x3, y1, y2, y3}]

Out[279]= {610., {x1 -> 80, x2 -> 150, x3 -> 0, y1 -> 1., y2 -> 1., 
  y3 -> 0.}}

In[278]:= NMaximize[{2 x1 + 3 x2 + 4 x3, 1.5 x1 + 3 x2 + 5 x3 <= 600, 
  280 x1 + 250 x2 + 400 x3 <= 60000, x1 <= 10^9*y1, x2 <= 10^9 y2, 
  x3 <= 10^9  y3, x1 >= 80 y1, x2 >= 80 y2, x3 >= 80 y3, x1 >= 0, 
  x2 >= 0, x3 >= 0, Element[{x1, x2, x3}, Integers], 
  y1 == 0 || y1 == 1, y2 == 0 || y2 == 1, y3 == 1 || y3 == 0}, {x1, 
  x2, x3, y1, y2, y3}]

During evaluation of In[278]:= NMaximize::bcons: The following constraints are not valid: {1.5 x1+3 x2+5 x3<=600,280 x1+250 x2+400 x3<=60000,x1<=1000000000 y1,x2<=1000000000 y2,x3<=1000000000 y3,x1>=80 y1,x2>=80 y2,x3>=80 y3,x1>=0,x2>=0,x3>=0,(x1|x2|x3)\[Element]Integers,y1==0||y1==1,y2==0||y2==1,y3==1||y3==0}. Constraints should be equalities, inequalities, or domain specifications involving the variables. >>

Out[278]= NMaximize[{2 x1 + 3 x2 + 4 x3, 1.5 x1 + 3 x2 + 5 x3 <= 600, 
  280 x1 + 250 x2 + 400 x3 <= 60000, x1 <= 1000000000 y1, 
  x2 <= 1000000000 y2, x3 <= 1000000000 y3, x1 >= 80 y1, x2 >= 80 y2, 
  x3 >= 80 y3, x1 >= 0, x2 >= 0, 
  x3 >= 0, (x1 | x2 | x3) \[Element] Integers, y1 == 0 || y1 == 1, 
  y2 == 0 || y2 == 1, y3 == 1 || y3 == 0}, {x1, x2, x3, y1, y2, y3}]

ref: https://reference.wolfram.com/language/guide/Optimization.html

解： 引入变量，当第人做工作时，取值为，否则为。则指派问题的数学模型为：

目标：

约束：

Python中，可以使用

scipy.optimize.linear_sum_assignment(cost_matrix)

https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.optimize.linear_sum_assignment.html#scipy.optimize.linear_sum_assignment

模型

用表示选修课程，表示不选。

目标:

约束:

1. 至少学习的课程数
2. 课程的先修课程要求。如：“数据结构”的选修课是“计算机编程”，表示若则必有，若，则对的值没有限制，可以用来表示。类似，“最优化方法”的先修课是“微积分”和“线性代数”，则有, ，表示为。

加上是变量的约束条件

Python

scipy.optimize

cvxpy (a Python-embedded modeling language for convex optimization problems)
https://stanford.edu/class/ee364a/lectures.html

Pyomo is a Python-based, open-source optimization modeling language with 
  a diverse set of optimization capabilities.

PuLP is an LP modeler written in python.
https://github.com/coin-or/pulp

import numpy as np

z = np.array([2, 3, 1])
a = np.array([[1, 4, 2], [3, 2, 0]])
b = np.array([8, 6])

x1_bound = x2_bound = x3_bound =(0, None)

from scipy import optimize

res = optimize.linprog(z, A_ub=-a, b_ub=-b,bounds=(x1_bound, x2_bound, x3_bound))

print(res)

#output:
#     fun: 7.0
# message: 'Optimization terminated successfully.'
#     nit: 2
#   slack: array([0., 0.])
#  status: 0
# success: True
#       x: array([0.8, 1.8, 0. ])

pulp在线性规划这一部分

https://www.jianshu.com/p/9be417cbfebb

import cvxpy
import numpy as np

# The data for the Knapsack problem
# P is total weight capacity of sack
# weights and utilities are also specified
P = 165
weights = np.array([23, 31, 29, 44, 53, 38, 63, 85, 89, 82])
utilities = np.array([92, 57, 49, 68, 60, 43, 67, 84, 87, 72])

# The variable we are solving for
selection = cvxpy.Bool(len(weights))

# The sum of the weights should be less than or equal to P
weight_constraint = weights * selection <= P

# Our total utility is the sum of the item utilities
total_utility = utilities * selection

# We tell cvxpy that we want to maximize total utility 
# subject to weight_constraint. All constraints in 
# cvxpy must be passed as a list
knapsack_problem = cvxpy.Problem(cvxpy.Maximize(total_utility), [weight_constraint])

# Solving the problem
knapsack_problem.solve(solver=cvxpy.GLPK_MI)

本节读完

17.