数学建模中非线性规划_线性规划数学建模案例

数学建模笔记（二）——非线性规划

一、基本介绍

线性规划与非线性规划的差距

线性规划是指模型中所有变量的次数只能是一次，非线性规划就是包含$x^a(a\ne 1),\sin (x),\cos (x),e^x,ln(x)$的模型。
$$\begin{gathered} minf(x)=x_1^2+x_2^2+x_3^2+8 \\ \{\begin{array}{l}{x}_1^2-x_2+x_3^2\ge 0\\ x_1+x_2^2+x_3^3\le 20\\ -x_1-x_2^2+2=0\\ x_2+2x_3^2=3\\ x_1,x_2,x_3\ge 0\end{array} \end{gathered}$$
在形式上和线性规划相当类似，解题步骤上也是类似，先确定决策变量，再写出目标函数和约束条件。

模型简介

虽然形式上与线性规划相当类似，但是求解却相当复杂。

• 因为线性规划有通用的求准确解的方法，一定可以求到最优解:happy:

• 非线性规划在数学上并没有求的精准数值解的通用方法

• 可也通过matlab的fmincon函数求的近似解

  非线性规划的解可以存在在可行域上的任意一点

二、求解方法

$fmincon$方法

$matlab$的$fmincon$函数可以求解：$[x,fval]=fmincon(fun,x_0,A,b,Aeq,beq,lb,ub,nonlcon)$

注意初始值$x_0$会影响求解结果，可以多次设不同的初始值进行比较，或者配合蒙特卡洛法进行求解。

$Scipy.optimize$模块

除了$matlab的fmincon$可以进行求解以外，我们通常还可以使用python的库$Scipy$的$Optimize$模块来进行求解。

$scipy.optimize$提供了多个用于非线性规划的方法：

• $brent()$:单变量无约束优化问题，用的是牛顿法/二分法。
• $fmin()$：多变量无约束优化问题，使用单纯形法，只需要利用函数值，不用函数的导数。
• $leatsq()$：非线性最小二乘拟合问题
• $minimize()$：约束优化问题，通过拉格朗日乘子法将约束问题转化成无约束问题

三、例题讲解

$matlab$例题

1.工地选址

某公司有6个建筑工地要开工，每个工地的位置（用平面坐标系a，b表示，距离单位：km） 及水泥日用量d（单位：t）由下表给出。目前有两个临时料场位于A（5，1），B（2，7）。 日储量各有==$20t$==。假设从料场到工地之间均有直线道路相连。

• 问题一：试制定每天的供应计划，即从A，B两料场分别向各工地运送多少吨水泥，使总的吨千米数最小？
• 问题二：为了进一步减少吨千米数，打算舍弃这两个临时料场，改建两个新的，日储量各为 $20t$。问应建在何处，节省的吨千米数为多大？

1、确定决策变量

• 第$i$个工地的坐标为$(a_i,b_i)$,水泥的日用量为$d_i,i=1,2,3\dots 6$;料场位置$(x_j,y_j)$,日储量为$e_j,j=1,2$,从料场$j$运往工地$i$的运送量为$X_{ij}$。
• 第一问中已确定料场的位置，所以距离已知质量，因此该题为线性规划问题，第二问则为非线性规划问题。
• 第一问决策变量为$X_{ij}$,第二问的决策变量为$x_j,y_j,X_{ij}$。

2、确定约束条件

每日的水泥用料量不能超过最大限额:
$$\sum _{i=1}^6{X}_{ij}\le e_j,j=1,2$$
每个工地的水泥量必须等于要求量:
$$\sum _{j=1}^2{X}_{ij}=d_i,i=1,2\dots 6$$
3、确定目标函数

要求总千米吨数最小，即运送量乘以运送里程的求和:
$$\sum _{i=1}^6\sum _{i=1}^2{X}_{ij}\sqrt{(x_j-a_i{)}^2+(y_j-b_i{)}^2},i=1,2\dots 6,j=1,2$$
连立，可得模型
$$\begin{gathered} minf=\sum _{i=1}^6\sum _{i=1}^2{X}_{ij}\sqrt{(x_j-a_i{)}^2+(y_j-b_i{)}^2} \\ s.t.\{\begin{array}{l}{\sum }_{i=1}^6{X}_{ij}\le e_j,j=1,2\\ {\sum }_{j=1}^2{X}_{ij}=d_i,i=1,2\dots 6\end{array} \end{gathered}$$

• 对于第一题，只有一个决策变量$X_{ij}$，目标函数中根号下的$x_j$和$y_j$是常数，目标函数和约束条件都是线性的，所以第一题是求解线性规划问题。
• 对于第二题，因为存在$\sqrt{(x_j-a_i{)}^2+(y_j-b_i{)}^2}$，所以是非线性的
  .

4,代码

clear;clc;
% 6个工地坐标 数据初始化
a=[1.25 8.75 0.5 5.75 3 7.25];
b=[1.25 0.75 4.75 5 6.5 7.75];
% 临时料场位置
x=[5 2];
y=[1 7];
% 6个工地水泥日用量
d=[3 5 4 7 6 11];
% 计算目标函数系数，即6工地与两个料场的距离，总共12个值
for i=1:6   % 对于6个工地
    for j=1:2       % 接收两个料场的供用
        l(i,j)=sqrt((x(j)-a(i))^2+(y(j)-b(i))^2);   % 距离
    end
end

%第一题的求解
f = [l(:,1);l(:,2)];    % 目标函数系数向量，总共12个值

% 不等式约束条件的变量系数和常数项
A = [1 1 1 1 1 1 0 0 0 0 0 0
     0 0 0 0 0 0 1 1 1 1 1 1];    
% 两个临时料场日储量
b = [20;20];

% 矩阵的行数是约束条件个数，列数是变量个数
% 等式约束的变量系数和常数项
Aeq = [eye(6),eye(6)];    % 两个单位矩阵横向拼成
beq=[d(1);d(2);d(3);d(4);d(5);d(6)];
% 所有变量下限全是0
Vlb=[0 0 0 0 0 0 0 0 0 0 0 0]; 
[x,fval]=linprog(f,A,b,Aeq,beq,Vlb);
x,fval

延续第一问设的$X_1$到$X_{12}$，代表料场到工地的运输水泥量；
第二问要求确定新料场的位置，所以要再添加4个变量，作为两个新料场的横纵坐标。
因此，第二问有16个变量，设定为$X_1$到$X_{16}$。
相应地，每个约束条件的系数矩阵中，也要有16个元素。

% 第二问中求新料场位置，所以两个料场的横纵坐标也是变量，所以多了4个变量
% 对新坐标没有不等式约束，所以其不等式约束条件里的系数为0
A2 = [1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0
      0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0];
B2 = [20;20];
% 对新坐标也没有等式约束，所以相应项也为0
Aeq2 = [eye(6),eye(6),zeros(6,4)];    % 两个单位矩阵和一个全0矩阵拼成
beq2=[3 5 4 7 6 11]';
Vlb2=[zeros(12,1);-inf;-inf;-inf;-inf];

$matlab的fmincon$必须赋初值,可以直接用第一问的结果作为初始值。

x0=[3 5 0 7 0 1 0 0 4 0 6 10 5 1 2 7]';  
[x2,fval2]=fmincon(@obj_f,x0,A2,B2,Aeq2,beq2,Vlb2)

本题只有目标函数有非线性函数，若约束条件中有非线性函数，可以使用nonlcon项，具体可以参考

https://ww2.mathworks.cn/help/optim/ug/fmincon.html?searchHighlight=fmincon&s_tid=srchtitle_fmincon_1#busog7r-nonlcon

若有条件，可使用蒙特卡罗法求一个近似解作为初始值
1、因为第2问模型中有16个变量，所以要给16个变量分别生成随机值，作为当前解；
2、判断这些当前解是否满足模型的约束条件
3、若满足，代入目标函数，求当前目标函数值
4、判断当前目标函数值是否比已求的较优函数值更好，若是，则替换掉较优函数值和对应的较优解
5、不断重复前5步，构成统计意义，求得较优解。

python例题

1,$brent$

例，求出函数$f(x)=x^2+8sin(2x+\pi )$的最大值。

scipy.optimize.brent(func, args=(), brack=None, tol=1.48e-08, full_output=0, maxiter=500)

$optimize.brent()$的主要参数

$optimize.brent()$的主要参数

代码如下

from scipy.optimize import brent, fmin_ncg, minimize
import numpy as np

# 1. Demo1：单变量无约束优化问题(Scipy.optimize.brent)
def objf(x):  # 目标函数
    fx = x**2 - 8*np.sin(2*x+np.pi)
    return fx

xIni = -5.0
xOpt= brent(objf, brack=(xIni,2))
print("xIni={:.4f}\tfxIni={:.4f}".format(xIni,objf(xIni))
print("xOpt={:.4f}\tfxOpt={:.4f}".format(xOpt,objf(xOpt)))

最后求得结果为

$$\begin{gathered} xIni=-5.0000,fxIni=29.3522 \\ xOpt=-0.7391,fxOpt=-7.4195 \end{gathered}$$

2,$fmin$

多变量无约束优化问题的算法很多，分类方式也很多。从使用者的角度来说可以分为：只使用目标函数值、使用导数（梯度下降法）、使用二阶导数。大体来说，使用导数的算法收敛较快，使用二阶导数收敛更快，但是收敛快也容易陷入局部最优。

$fmin()$ 函数是 SciPy.optimize 模块中求解多变量无约束优化问题（最小值）的首选方法，采用下山单纯性方法。下山单纯性方法又称 Nelder-Mead 法，只使用目标函数值，不需要导数或二阶导数值，是最重要的多维无约束优化问题数值方法之一。

scipy.optimize.fmin(func, x0, args=(), xtol=0.0001, ftol=0.0001, maxiter=None, maxfun=None, full_output=0, disp=1, retall=0, callback=None, initial_simplex=None)

案例：

这题是一个经典函数，用来测试优化算法性能的非凸函数。叫$Rosenbrock$山谷函数
$$f(x,y)=(1-x{)}^2+100(y-x^2{)}^2$$
代码中通过一个二元列表代替了x和y两个变量。

from scipy.optimize import brent, fmin, minimize
import numpy as np

# 2. Demo2：多变量无约束优化问题(Scipy.optimize.brent)
# Rosenbrock 测试函数
def objf2(x):  # Rosenbrock benchmark function
    fx = sum(100.0 * (x[1:] - x[:-1] ** 2.0) ** 2.0 + (1 - x[:-1]) ** 2.0)
    return fx

xIni = np.array([-2, -2])
xOpt = fmin(objf2, xIni)
print("xIni={:.4f},{:.4f}\tfxIni={:.4f}".format(xIni[0],xIni[1],objf2(xIni)))
print("xOpt={:.4f},{:.4f}\tfxOpt={:.4f}".format(xOpt[0],xOpt[1],objf2(xOpt)))

今天的文章数学建模中非线性规划_线性规划数学建模案例分享到此就结束了，感谢您的阅读。