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)
? 线性规划中对于初始值 x0 的选取? 常重要 , 因为? 线性规划的算法求解出来的 是? 个局部优解 . (线性规划不存在这个问题 )
如果要求全局优解! 有两种思路 :
给定不同的初始值 , 在?? 找到? 个优解 ;
先? 蒙特卡罗模拟, 得到? 个蒙特卡罗解 , 然后将这个解作为初始值来求优解
"option" 选项可以给定求解的算法, ? 共有四种 :
- interior.point (内点法)
- sqp(序列? 次规划法)
- active.set (有效集法)
- trust-region.reflective (信赖域 反射算法) .
不同的算法有其各? 的优缺点和适? 情况 , 我们可以改变求解的算法来看求解的结果是否变好了 . 如何改变求解的算法请参考代码演示 . (数模? 赛中可以都尝试 下, 这可以体现出稳健性, 也是你的优点)
"@fun" 表示? 标函数
我们要编写? 个独? 的 m ? 件来存? 标函数 :
- function f = fun ( x)
- f=...
- end
注 fun 可以任意取名 , 只 要符合 Matlab 命名规范 , 保存的 m ? 件 也是这个名 .
f 也可以任意取名, 但返回的 f 和 函数内部的 f 得 完全? 致;
这? 的 x 实际上是表示决策变量的向量 , 其行列方向取决于初始值.
调? 函数 : fmincon(@fun,...) 求解 .
"@nonlfun" 表示? 线性部分的约束 , 我们同样得编写? 个独? 的 m ? 件储存非线性约束条件 :
nonlfun 同样可任意取名 , 别和上? 的 fun 相同即可, 保存的 m 文件也得是 这个 名字.
c,ceq 中可能有多个约束 因此我们写成列向量的形式 ;
若不存在? 线性不等式约束 , 则可以令 C = []
调? 函数 fmincon(.......,@nonlfun.options) 求解 .
注意要把下标改写为扩号 , 例如 \(f=x_{1}^2+3x_{2}\)名 写成 Matlab 能识别的就应该为 : \(f=x(1)^2+3*x(2)\)
若 不存在某种约束 , 则 可? "[]" 替代 , 若后? 全为 "[ ]" 且 不指定 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
典型例题
来源: http://www.bubuko.com/infodetail-3438396.html