MATLAB使用蒙特卡洛算法实例求非线性规划

蒙特卡洛方法是一种利用计算机的随机数理论模拟实际的情况的一种方法。今天主要是以实例讲解蒙特卡洛方法的MATLAB编程实现求解非线性规划和非线性整数规划。

实例1

首先使用fmincon函数求解非线性规划，fmincon函数常用语句的语法如下：

[x,fval]=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
x的返回值是决策向量x的取值，fval的返回值是目标函数f(x)的取值
fun是用M文件定义的函数f(x),代表了(非)线性目标函数
x0是x的初始值
A,b,Aeq,beq定义了线性约束 ,如果没有线性约束，则A=[],b=[],Aeq=[],
beq=[]
lb和ub是变量x的下界和上界，如果下界和上界没有约束，则lb=[],ub=[],
也可以写成lb的各分量都为 -inf,ub的各分量都为inf
nonlcon是用M文件定义的非线性向量函数约束
options定义了优化参数，不填写表示使用Matlab默认的参数设置

主程序

clc;
clear all;
[x,y]=fmincon('fun1',rand(3,1),[],[],[],[],zeros(3,1),[],'fun2')

fun1.m自定义函数

function f=fun1(x);


f=x(1).^2+x(2).^2+x(3).^2+8;


end

fun2.m自定义函数

function [g,h]=fun2(x);
g=[-x(1).^2+x(2)-x(3).^2;
    x(1)+x(2).^2+x(3).^3-20];
h=[-x(1)-x(2).^2+2;
    x(2)+2*x(3).^2-3];


end

运行结果

x =
    0.5522
    1.2033
    0.9478
y =


   10.6511

蒙特卡洛求解带等式约束的非线性规划程序

主程序

clc;
clear all;
rand('state',sum(clock));%初始化随机数发生器
f0=inf;
x0 = [];
num = 1e6;
tic%计时开始
for i=1:num
    x3=0 + 20*rand(1,1);%随机产生初始解
    %% 尝试3 ：将等式约束转化在随机生成的过程中
    % h=[-x(1)-x(2).^2+2;
    %     x(2)+2*x(3).^2-3];
    x2  = -2*x3^2+3;
    x1 = -x2^2+2;
    x = [x1 x2 x3];
    [f,g]=mengte6(x);%调用自定义函数计算
    
    if (sum(g<=0)==2)
        if f0>=f  %求最小值 如果当前值更优，则更新值
            x0=x;
            f0=f;
           
        end
    end
end
toc%计时结束
fprintf('min f(x) 在x1 = %f x2 = %f  x3 = %f处取得最小值：%f\n',x0(1),x0(2),x0(3),f0);

mengte6.m自定义函数

function [f,g,h]=mengte6(x)
%% f是目标函数  g(x)<=0
%尝试1 h(x)==0  
%尝试2 h(x)>=0 h(x)<=0 
f=x(1).^2+x(2).^2+x(3).^2+8;
g=[-x(1).^2+x(2)-x(3).^2;
    x(1)+x(2).^2+x(3).^3-20;
  ];
% h=[-x(1)-x(2).^2+2;
%     x(2)+2*x(3).^2-3];
end

运行结果

历时 1.154699 秒。
min f(x) 在x1 = 0.552347 x2 = 1.203185  x3 = 0.947844处取得最小值：10.651148

实例2

首先使用fmincon函数求解非线性规划

主程序

clc;
clear all;
close all;
A = [1 1 1 1 1;1 2 2 1 6;2 1 6 0 0;0 0 1 1 5];
b = [400;800;200;200];
Aeq = [];
beq = [];
lb = zeros(5,1);
ub = zeros(5,1)+99;
x0 = rand(5,1);
[x,y]=fmincon(@fun11,x0,A,b,Aeq,beq,lb,ub,@fun22);
x = x
y_max = -y

自定义函数fun11.m

function y = fun11(x)
x1 = x(1);
x2 = x(2);
x3 = x(3);
x4 = x(4);
x5 = x(5);
y = x1^2+x2^2+3*x3^2+4*x4^2+2*x5^2-8*x1-2*x2-3*x3-x4-2*x5;
y = -y;
end

自定义函数fun22.m

function [g,h] = fun22(x)
%% g(x)<=0 h(x)=0  非线性约束条件
g = [];
h = [];
end

运行结果

Local minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.
<stopping criteria details>
x =
    0.0000
    0.0000
   33.3333
   99.0000
   13.5333
y_max =
   4.2678e+04

蒙特卡洛求解非线性规划程序

主程序

clc;
clear all;
rand('state',sum(clock));%初始化随机数发生器
f0=-inf;
x0 = [];
num = 1e6;
tic%计时开始
for i=1:num
    x=0 + 99*rand(1,5);%随机产生初始解
    [f,g]=mengte7(x);%调用自定义函数计算
    if (sum(g<=0)==4)
        if f0<=f  %求最大值 如果当前值更优，则更新值
            x0=x;
            f0=f;
           
        end
    end
end
toc%计时结束
fprintf('max f(x) 在x1 = %f x2 = %f  x3 = %f x4 = %f  x5 = %f处取得最大值：%f\n',x0(1),x0(2),x0(3),x0(4),x0(5),f0);

自定义函数mengte7.m

function [f,g]=mengte7(x)
%% f是目标函数  g(x)<=0
x1 = x(1);
x2 = x(2);
x3 = x(3);
x4 = x(4);
x5 = x(5);
f = x1^2+x2^2+3*x3^2+4*x4^2+2*x5^2-8*x1-2*x2-3*x3-x4-2*x5;


g=[x1+x2+x3+x4+x5-400;
   x1+2*x2+2*x3+x4+6*x5-800;
   2*x1+x2+6*x3-200;
   x3+x4+5*x5-200;
  ];


end

运行结果

历时 1.021252 秒。
max f(x) 在x1 = 45.809604 x2 = 90.187886  x3 = 0.858626 x4 = 98.697572  x5 = 2.125203处取得最大值：48556.085563
>> 

实例3

蒙特卡洛求解非线性规划程序

主程序

clc;
clear all;
rand('state',sum(clock));%初始化随机数发生器
f0=-inf;
x0 = [];
num = 1e7;
tic%计时开始
for i=1:num
%     x=0 + 99*rand(1,5);%随机产生初始解
      x = randi([0 99],1,5);
    [f,g]=mengte7(x);%调用自定义函数计算
    if (sum(g<=0)==4)
        if f0<=f  %求最大值 如果当前值更优，则更新值
            x0=x;
            f0=f;
           
        end
    end
end
toc%计时结束
fprintf('max f(x) 在x1 = %f x2 = %f  x3 = %f x4 = %f  x5 = %f处取得最大值：%f\n',x0(1),x0(2),x0(3),x0(4),x0(5),f0);

自定义函数mengte7.m

function [f,g]=mengte7(x)
%% f是目标函数  g(x)<=0
x1 = x(1);
x2 = x(2);
x3 = x(3);
x4 = x(4);
x5 = x(5);
f = x1^2+x2^2+3*x3^2+4*x4^2+2*x5^2-8*x1-2*x2-3*x3-x4-2*x5;


g=[x1+x2+x3+x4+x5-400;
   x1+2*x2+2*x3+x4+6*x5-800;
   2*x1+x2+6*x3-200;
   x3+x4+5*x5-200;
  ];


end

运行结果

历时 14.614175 秒。
max f(x) 在x1 = 50.000000 x2 = 96.000000  x3 = 0.000000 x4 = 99.000000  x5 = 17.000000处取得最大值：50773.000000
>> 

本文内容来源于网络，仅供参考学习，如内容、图片有任何版权问题，请联系处理，24小时内删除。

作 者 | 郭志龙
编 辑 | 郭志龙
校 对 | 郭志龙