数学规划模型——非线性规划问题

Matlab 中⾮线性规划的标准型

\[min \ f(x) \\ \ \\ s.t. \ \begin{cases} & \text{} AX<=b,Aeg*x=beg \quad 线性约束\\ & \text{} c(x)<=0,ceg(x)=0\quad 非线性约束 \\ & \text{}lb<=x<=ub \quad 上下界约束(也可以当成不等式约束) \end{cases} \]

练习题

Matlab求解⾮线性规划的函数

[x,fval] = fmincon(@fun,x0,A,b,Aeq,beq,lb,ub,@nonlfun,option)

  1. ⾮线性规划中对于初始值 x0 的选取⾮常重要 , 因为⾮线性规划的算法求解出来的 是⼀个局部优解 。 (线性规划不存在这个问题 )

  2. 如果要求全局优解! 有两种思路 :

    • 给定不同的初始值 , 在⾥⾯找到⼀个优解 ;
    • 先⽤蒙特卡罗模拟, 得到⼀个蒙特卡罗解 , 然后将这个解作为初始值来求优解

  3. "option" 选项可以给定求解的算法, ⼀共有四种 :

    • interior.point (内点法)

    • sqp(序列⼆次规划法)

    • active.set (有效集法)

    • trust-region.reflective (信赖域 反射算法) 。

  4. 不同的算法有其各⾃的优缺点和适⽤情况 , 我们可以改变求解的算法来看求解的结果是否变好了 。 如何改变求解的算法请参考代码演示 。 (数模⽐赛中可以都尝试 下, 这可以体现出稳健性, 也是你的优点)

  5. "@fun" 表示⽬标函数
    我们要编写⼀个独⽴的 m ⽂件来存⽬标函数 :

    function f = fun ( x) 
    f=...
end
    • 注 fun 可以任意取名 , 只 要符合 Matlab 命名规范 , 保存的m ⽂件 也是这个名 。
    • f 也可以任意取名, 但返回的 f 和 函数内部的f得 完全⼀致;
    • 这⾥的x实际上是表示决策变量的向量 ,其行列方向取决于初始值。
    • 调⽤函数 : fmincon(@fun,...) 求解 。

  6. " @nonlfun " 表示⾮线性部分的约束 , 我们同样得编写⼀个独⽴的 m ⽂件储存非线性约束条件 :

    • nonlfun 同样可任意取名 , 别和上⾯的fun相同即可, 保存的m文件也得是 这个 名字。
    • c,ceq中可能有多个约束 因此我们写成列向量的形式 ;
    • 若不存在⾮线性不等式约束 , 则可以令 C = []
    • 调⽤函数 fmincon(.......,@nonlfun.options) 求解 。

  7. 注意要把下标改写为扩号 , 例如\(f=x_{1}^2+3x_{2}\)名 写成 Matlab 能识别的就应该为 : \(f=x(1)^2+3*x(2)\)

  8. 若 不存在某种约束 , 则 可⽤ " [] "替代 , 若后⾯全为 "[ ] " 且 不指定option(使⽤默认的求解⽅法) , 则"[]"也可以省略掉 。

练习题代码

code.m

% 非线性规划的函数
% [x,fval] = fmincon(@fun,x0,A,b,Aeq,beq,lb,ub,@nonlfun,option)
% x0表示给定的初始值（用行向量或者列向量表示），必须得写
% A b表示线性不等式约束
% Aeq beq 表示线性等式约束
% lb ub 表示上下界约束
% @fun表示目标函数
% @nonlfun表示非线性约束的函数
% option 表示求解非线性规划使用的方法
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文件和nonfun1.m文件都必须在当前文件夹目录下
fval = -fval
% 一个值得讨论的地方，能不能把线性不等式约束Ax <= b也写到nonlfun1函数中？
% 先把nonlfun1中的c改为下面这样：
% c = [(x(1)-1)^2-x(2); 
%        -2*x(1)+3*x(2)-6];
%  [x,fval] = fmincon(@fun1,x0,[],[],[],[],[],[],@nonlfun1)
% 结果也是可以计算出来的，但并不推荐这样做~

%% 使用其他算法对例题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文件和nonfun1.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*x1(i)+3*x2(i)-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  



%% 例题二的求解
x0 = [1 1 1];  %任意给定一个初始值 
lb = [0 0 0];  % 决策变量的下界
[x,fval] = fmincon(@fun2,x0,[],[],[],[],lb,[],@nonlfun2)  % 注意 fun2.m文件和nonfun2.m文件都必须在当前文件夹目录下
% x =
%          0.552167405729277          1.20325915507969         0.947824046150443
% fval =
%           10.6510918606939



%% 使用蒙特卡罗的方法来找初始值(推荐）
clc,clear;
n=1000000; %生成的随机数组数
x1= unifrnd(0,2,n,1);   % 生成在[0,2]之间均匀分布的随机数组成的n行1列的向量构成x1
x2 = sqrt(2-x1);  % 根据非线性等式约束用x1计算出x2
x3 = sqrt((3-x2)/2); % 根据非线性等式约束用x2计算出x3
fmin=+inf; % 初始化函数f的最小值为正无穷（后续只要找到一个比它小的我们就对其更新）
for i=1:n
    x = [x1(i), x2(i), x3(i)];  %构造x向量, 这里千万别写成了：x =[x1, x2, x3]
    if (-x(1)^2+x(2)-x(3)^2<=0) & (x(1)+x(2)^2+x(3)^2-20<=0)   % 判断是否满足条件
        result =sum(x.*x) + 8 ;  % 如果满足条件就计算函数值
        if  result  < fmin  % 如果这个函数值小于我们之前计算出来的最小值
            fmin = result;  % 那么就更新这个函数值为新的最小值
            x0 = x;  % 并且将此时的x1 x2 x3更新为初始值
        end
    end
end
disp('蒙特卡罗选取的初始值为：'); disp(x0)
lb = [0 0 0];  % 决策变量的下界
[x,fval] = fmincon(@fun2,x0,[],[],[],[],lb,[],@nonlfun2)  % 注意 fun2.m文件和nonfun2.m文件都必须在当前文件夹目录下

%% 例题三的求解(蒙特卡罗模拟那一讲的例题)
clear;clc
% 蒙特卡罗模拟得到的最大值为3445.6014
% 最大值处x1 x2 x3的取值为：
%           22.5823101903968          12.5823101903968          12.1265223966757
A = [1 -2 -2;  1 2 2];  b = [0 72];
x0 = [ 22.58   12.58  12.13];
Aeq = [1 -1 0]; beq = 10;
lb = [-inf 10 -inf];  ub = [inf 20 inf];  
[x,fval] = fmincon(@fun3,x0,A,b,Aeq,beq,lb,ub,[])  % 注意没有非线性约束，所以这里可以用[]替代，或者干脆不写
fval = -fval

fun1.m

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

fun2.m

function f = fun2(x)
    %     f = x(1)^2+x(2)^2 +x(3)^2+8 ; 
    f = sum(x.*x) + 8;  % 可别忘了x实际上是一个向量，我们可以使用矩阵的运算符号对其计算
end

fun3.m

function f = fun3(x)
    f = -prod(x);  % 可别忘了x实际上是一个向量（prod表示连乘符号，用法和sum类似）
end

nonlfun1.m

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

nonlfun2.m

function [c,ceq] = nonlfun2(x)
    % 非线性不等式约束
    c = [-x(1)^2+x(2)-x(3)^2;   % 一定要注意写法的规范，再次强调这里的x是一个向量！不能把x(1)写成x1
            x(1)+x(2)^2+x(3)^2-20];
    % 非线性等式约束
    ceq = [-x(1)-x(2)^2+2;
                x(2)+2*x(3)^2-3]; 
end

典型例题

posted @ 2020-02-27 19:41  pxlsdz  阅读(2717)  评论(0)  编辑  收藏  举报