《零基础数学建模》——非线性规划

​ 前言

本文大部分是对于数学建模清风老师的课程学习总结归纳而来，我的理解可能有错误，大家发现错误可以在评论区批评指正，课程地址：《数学建模清风》
  本文的前置文章《零基础数学建模》——线性规划

---

一、非线性规划模型的概述

非线性规划（Nonlinear Programming）
  当⽬标函数和或者约束条件中有⼀个是决策变量 $x$ 的非线性表达式 , 那么此时的数学规划问题就属于非线性规划 。
  解决非线性规划要比线性规划困难得多 , ⽬前没有通⽤算法,⼤多数算法都是在选定决策变量的初始值后 ,通过一定的搜索⽅法寻求最优的决策变量 。

二、$Matlab$标准型

1. 非线性规划的标准型

$$\begin{gathered} minf(x) \\ s.t.\{\begin{array}{c}Ax\le b,Aeq\cdot x=beq(linearconstraints)\\ c(x)\le 0,ceq(x)=0（nonlinearconstraints）\\ lb\le x\le ub\end{array} \end{gathered}$$
PS：可能只对部分决策变量有约束。

2. 非线性规划的函数

$$[x,fval]=fmincon(@fun,x_0,A,b,Aeq,beq,lb,ub,@nonlfun,option)$$
注意：

1. 非线性规划中对于初始值$x_0$的选取非常重要，因为非线性规划的算法求解出来的是一个局部最优解。（线性规划不存在这个问题）
2. 如果要求“全局最优解”，有两个思路：给定不同的初始值，在里面找到一个最优解；先用蒙特卡罗模拟，得到一个蒙特卡罗解，然后将这个解作为初始值来求最优解。
3. "option"选项可以给定求解的算法，一共有四种：$interior-point$（内点法）、$sqp$（序列二次规划法）、$active-set$（有效集法）以及$trust-region-reflctive$（信赖与反射算法）。
4. 不同的算法有其各自的优缺点和使用情况，我们可以改变求解的算法来看求解的结果是否变好了。如何改变求解的算法请参考下文的经典例题代码部分。（比赛中可以都尝试下，这可以体现出稳健性，也是团队的优点）
5. "@fun"表示目标函数，我们要编写一个独立的“m文件”储存目标函数：

6. “@nonlfun”表示非线性部分的约束，我们同样需要编写一个独立的“m文件”储存非线性约束条件：

7. 注意要把下标改写为括号，例如：$f=x_1^2+3x_2$写成Matlab能识别的形式应该为：$f=x(1{)}^2+3\ast x(2)$。
8. 若不存在某种约束，则可用"[]“替代，若后面全为”[]“且不指定"option”（使用默认的求解方法），则"[]"可以省略掉。

三、经典例题

1. 例题1的求解

主函数求解：

clear;clc
format long g   %可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）

%% 例题1的求解
% max f(x) = x1^2 +x2^2 -x1*x2 -2x1 -5x2
% s.t. -(x1-1)^2 +x2 >= 0 ;  2x1-3x2+6 >= 0
x0 = [0 0];  %任意给定一个初始值 
A = [-2 3]; b = 6;
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1)  % 注意 fun1.m文件和nonlfun1.m文件都必须在当前文件夹目录下
fval = -fval

目标函数

function f = fun1(x)
    % 注意：这里的f实际上就是目标函数，函数的返回值也是f
    % 输入值x实际上就是决策变量，由x1和x2组成的向量
    % fun1是函数名称，到时候会被fmincon函数调用, 可以任意取名
    % 保存的m文件和函数名称得一致，也要为fun1.m
%      max  f(x) = x1^2 +x2^2 -x1*x2 -2x1 -5x2
    f = -x(1)^2-x(2)^2 +x(1)*x(2)+2*x(1)+5*x(2) ; 
end

非线性约束

function [c,ceq] = nonlfun1(x)
    % 注意：这里的c实际上就是非线性不等式约束，ceq实际上就是非线性等式约束
    % 输入值x实际上就是决策变量，由x1和x2组成的一个向量
    % 返回值有两个，一个是非线性不等式约束c，一个是非线性等式约束ceq
    % nonlfun1是函数名称，到时候会被fmincon函数调用, 可以任意取名，但不能和目标函数fun1重名
    % 保存的m文件和函数名称得一致，也要为nonlfun1.m
%     -(x1-1)^2 +x2 >= 0 
   c = [(x(1)-1)^2-x(2)];   % 千万別写成了: (x1-1)^2 -x2
   ceq = [];  % 不存在非线性等式约束，所以用[]表示
end

主函数更换求解算法以及选取初始值方法

%% 使用其他算法对例题1求解
% edit fmincon  % 查看fmincon的“源代码”
% Matlab2017a默认使用的算法是'interior-point' 内点法
% 使用interior point算法 （内点法）
option = optimoptions('fmincon','Algorithm','interior-point')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  
fval = -fval
% 使用SQP算法 （序列二次规划法）
option = optimoptions('fmincon','Algorithm','sqp')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  
fval = -fval   %得到-4.358，远远大于内点法得到的-1,猜想是初始值的影响
% 改变初始值试试
x0 = [1 1];  %任意给定一个初始值 
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  % 最小值为-1，和内点法相同（这说明内点法的适应性要好）
fval = -fval  
% 使用active set算法 （有效集法）
option = optimoptions('fmincon','Algorithm','active-set')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)
fval = -fval  
% 使用trust region reflective (信赖域反射算法)
option = optimoptions('fmincon','Algorithm','trust-region-reflective')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  
fval = -fval  
% this algorithm does not solve problems with the constraints you have specified. 
% 这说明这个算法不适用我们这个约束条件，所以以后遇到了不能求解的情况，记得更换其他算法试试！！！

%% 选取初始值得到的结果可能会不满足限定条件，出现了一个Bug 因此选择的初始值很重要
x0 = [40.8, 10.8];
option = optimoptions('fmincon','Algorithm','interior-point')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  
fval = -fval  
% https://cn.mathworks.com/help/optim/ug/fmincon.html

%% 生成不同的随机初始值来优化代码，有一定几率会触发上面那个Bug，因此不推荐
n = 10;  % 重复n次
Fval = +inf; X = [0,0];  %初始化最优的结果
A = [-2 3]; b = 6;
for i = 1:n
    x0 = [rand()*10 , rand()*10];  %用随机数生成一个初始值（随机数的范围自己根据题目条件设置） 
    [x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option); % 注意 fun1.m文件和nonlfun1.m文件都必须在当前文件夹目录下
    if fval < Fval  % 如果找到了更小的值，那么就代替最优的结果
        Fval = fval;
        X = x;
    end
end
Fval = -Fval
X

%% 使用蒙特卡罗的方法来找初始值(推荐）
clc,clear;
n=10000000; %生成的随机数组数
x1=unifrnd(-100,100,n,1);  % 生成在[-100,100]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=unifrnd(-100,100,n,1);  % 生成在[-100,100]之间均匀分布的随机数组成的n行1列的向量构成x2
fmin=+inf; % 初始化函数f的最小值为正无穷（后续只要找到一个比它小的我们就对其更新）
for i=1:n
    x = [x1(i), x2(i)];  %构造x向量, 这里千万别写成了：x =[x1, x2]
    if ((x(1)-1)^2-x(2)<=0)  & (-2*x(1)+3*x(2)-6 <= 0)     % 判断是否满足条件
        result = -x(1)^2-x(2)^2 +x(1)*x(2)+2*x(1)+5*x(2) ;  % 如果满足条件就计算函数值
        if  result  < fmin  % 如果这个函数值小于我们之前计算出来的最小值
            fmin = result;  % 那么就更新这个函数值为新的最小值
            x0 = x;  % 并且将此时的x1 x2更新为初始值
        end
    end
end
disp('蒙特卡罗选取的初始值为：'); disp(x0)
A = [-2 3]; b = 6;
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1)
fval = -fval  

2. 选址问题

  第一问是线性规划问题，我们这里只关注第二问的求解。

  记工地$i$的坐标为$(a_i,b_i),i=1,2,\cdots ,6$，水泥日用量为$d_i$；
  记新料场$j$的坐标为$(x_j,y_j),j=1,2$，日储量为20;
  从料场$j$往工地$i$运送的水泥为$x_{ij}$吨，则
$$\begin{gathered} minz=\sum _{i=1}^6\sum _{j=1}^2{x}_{ij}\ast \sqrt{(a_i-x_j{)}^2+(b_i-y_j{)}^2} \\ s.t.\{\begin{array}{c}{\sum }_{j=1}^2{x}_{ij}=d_i,i=1,2,\cdots ,6\\ {\sum }_{i=1}^6{x}_{ij}\le 20,j=1,2\\ x_{ij}\ge 0,i=1,2,\cdots ,6,j=1,2\end{array} \end{gathered}$$
转换为$matlab$标准型：

代码结果：

%% 选址问题
clear;clc
format long g   %可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）
% % (1) 系数向量（原来线性规划问题的写法,我们只需要在此基础上改动一点就可以了）
% a=[1.25  8.75  0.5  5.75  3  7.25];  % 工地的横坐标
% b=[1.25  0.75  4.75	5  6.5  7.25];   % 工地的纵坐标
% x = [5  2];  % 料场的横坐标
% y = [1  7];  % 料场的纵坐标
% c = [];  % 初始化用来保存工地和料场距离的向量 (这个向量就是我们的系数向量）
% for  j =1:2
%     for i = 1:6
%         c = [c;  sqrt( (a(i)-x(j))^2 + (b(i)-y(j))^2)];  % 每循环一次就在c的末尾插入新的元素
%     end
% end
% (2) 不等式约束
A =zeros(2,16);  % 注意这里要改成16
A(1,1:6) = 1;
A(2,7:12) = 1;
b = [20,20]';
% (3) 等式约束
Aeq = zeros(6,16);  % 注意这里要改成16
for i = 1:6
    Aeq(i,i) = 1;  Aeq(i,i+6) = 1;
end
beq = [3 5 4 7 6 11]';  % 每个工地的日需求量
%（4）上下界
lb = zeros(16,1);
% lb = [zeros(12,1); -inf*ones(4,1)];  两个新料场坐标的下界可以设为-inf

% 进行求解
% 注意哦，这里我们只尝试了这一个初始值，大家可以试试其他的初始值，有可能能够找到更好的解。
x0 = [3 5 0 7 0 1 0 0 4 0 6 10 5 1 2 7];  % 用第一问的结果作为初始值
[x,fval] = fmincon(@fun5,x0,A,b,Aeq,beq,lb)  % 注意没有非线性约束，所以这里可以用[]替代，或者干脆不写
reshape(x(1:12),6,2)  % 将x的前12个元素变为6行2列便于观察（reshape函数是按照列的顺序进行转换的，也就是第一列读完，读第二列，即x1对应x_1,1，x2对应x_2,1）
% 新坐标（5.74，4.99） （7.25，7.25）
% fval =
%           89.9231692432933
% 第一问的fval =
%           135.281541790676
135.281541790676 - 89.9231692432933  %  45.3583725473827

function f = fun5(xx)  % 注意为了避免和下面的x同号，我们把决策变量的向量符号用xx表示（注意xx的长度为16）
    a=[1.25  8.75  0.5  5.75  3  7.25];  % 工地的横坐标
    b=[1.25  0.75  4.75	5  6.5  7.25];   % 工地的纵坐标
    x = [xx(13)  xx(15)];  % 新料场的横坐标
    y = [xx(14)  xx(16)];  % 新料场的纵坐标
    c = [];  % 初始化用来保存工地和料场距离的向量 (这个向量就是我们的系数向量）
    for  j =1:2
        for i = 1:6
            c = [c;  sqrt( (a(i)-x(j))^2 + (b(i)-y(j))^2)];  % 每循环一次就在c的末尾插入新的元素
        end
    end
    % 下面我们要求吨千米数，注意c是列向量，我们计算非线性规划时给定的初始值x0是行向量
    f = xx(1:12) * c;
end