非线性规划matlab实现

非线性规划求得的是局部最优解，若要求“全局最优解”，可先用蒙特卡洛模拟，求得蒙特卡洛解作为初始值，再用非线性规划。

%% 非线性规划的函数
% [x,fval] = fmincon(@fun,x0,A,b,Aeq,beq,lb,ub,@nonlfun,option)
% @fun表示目标函数，写到m文件中，fun(x)传入的x是决策变量的向量
% x0表示给定的初始值（用行向量或者列向量表示），必须得写,推荐蒙特卡洛法
% 蒙特卡洛思想：随机生成n组x1和x2，然后依次计算每组x下的f(x),并使用更优的x进行更新迭代，n足够大时，以求得一个较优的解
n=10000000; %生成的随机数组数
% A b表示线性不等式约束（仅有一次项）
% Aeq beq 表示线性等式约束（仅有一次项）
% lb ub 表示上下界约束
% @nonlfun表示非线性约束的函数（有非一次项），写到m文件中，有两个返回值C和Ceq，分别为非线性不等式约束和非线性等式约束。不存在时为[]
% 非线性不等式约束c = [-x(1)^2+x(2)-x(3)^2;x(1)+x(2)^2+x(3)^2-20];
% 非线性等式约束ceq = [-x(1)-x(2)^2+2;x(2)+2*x(3)^2-3]; 
% 一定要注意写法的规范，再次强调这里的x是一个向量！不能把x(1)写成x1
% option 表示求解非线性规划使用的方法，一共4种：interior point算法 （内点法）、SQP算法 （序列二次规划法）、active set算法 （有效集法）、trust region reflective (信赖域反射算法)


%% 例题1的求解，两个决策变量x(1),x(2)
% max f(x) = x1^2 +x2^2 -x1*x2 -2x1 -5x2
% s.t. -(x1-1)^2 +x2 >= 0 非线性不等式 ;  2x1-3x2+6 >= 0 线性不等式



%% 使用蒙特卡罗的方法来找初始值x0(推荐）
clc,clear;
% 思想：随机生成n组x1和x2，然后依次计算每组x下的f(x),并使用更优的x进行更新迭代，n足够大时，以求得一个较优的解
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) ;  % 如果满足非线性不等式约束条件,就计算函数值f(x)
        if  result  < fmin  % 如果这个函数值f(x)小于我们之前计算出来的最小值fmin
            fmin = result;  % 那么就更新这个函数值为新的最小值
            x0 = x;  % 并且将此时的x1 x2更新为初始值
        end
    end
end
disp('蒙特卡罗选取的初始值为：'); disp(x0)
A = [-2 3]; b = 6;




%% 使用4种算法option分别对例题1求解
% 比赛时可以把所有算法都试一下，进行对比实验

% 使用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. 
% 这说明这个算法不适用我们这个约束条件，所以以后遇到了不能求解的情况，记得更换其他算法试试！！！

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 转min形加负号，注意x是向量用x(1)而不是x1
    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