首页 > 其他 > 详细

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

时间:2020-02-27 20:09:28      阅读:95      评论:0      收藏:0      [点我收藏+]

title: 数学规划模型——非线性规划问题
date: 2020-02-26 20:10:37
categories: 数学建模
tags: [MATLAB, 数学规划模型]


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

典型例题

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

原文:https://www.cnblogs.com/pxlsdz/p/12373973.html

(0)
(0)
   
举报
评论 一句话评论(0
关于我们 - 联系我们 - 留言反馈 - 联系我们:wmxa8@hotmail.com
© 2014 bubuko.com 版权所有
打开技术之扣,分享程序人生!