数学建模更新12(非线性规划)
非线性规划
- 一.非线性规划的标准型
- 1.概述
- 2.例题
- 3.MATLAB的非线性规划的函数
- 4.基础例题
- 【1】默认方法求例1
- 【2】用其他方法求例2
- 【3】改变初始值的影响
- 【4】使用蒙特卡罗的方法来找初始值(推荐)
- 【5】其他例题的求解
- 二.例题
- 1.选址问题
- 2.飞行管理问题
一.非线性规划的标准型
1.概述
2.例题
3.MATLAB的非线性规划的函数
[x,fval] = fmincon(@fun,x0,A,b,Aeq,beq,lb,ub,@nonlfun,option)
解释:
- 非线性规划中对于初始值xxx0的选取非常重要,因为非线性规划的算法求解出来的是一个局部最优解
- 如果要求“全局最优解”有两种思路:①给定不同的初始值,在里面得到一个最优解,②先用蒙特卡洛模拟,得到一个蒙特卡洛解,然后将这个解作为初始值求最优解
- "option""option""option"选项可以给定求解的算法,一共有四种:interior−pointinterior-pointinterior−point(内点法),sqpsqpsqp(序列二次规划),active−setactive-setactive−set(有效算法)以及trust−region−reflectivetrust-region-reflectivetrust−region−reflective(信赖域反射算法)
- 不同算法有其各自的优缺点和适用情况,我们可以改变求解的算法来看求解的结果是否变好了,(数模比赛中可以尝试下,这可以体现出稳健性)
- “@fun”“@fun”“@fun”表示目标函数,我们要编写一个独立的“m文件”储存目标函数
function F = fun(x)f=
end
①fun可以任意取值,符合MATLAB的命名规范就行,注意保存的m文件就是这个名字
②f也可以任意命名,但是返回的f和函数内部的f一致
③这里的x实际上是表示决策变量的向量,其行列方向取决于初始值x
④调用函数:fmincon(@fun,……)
求解
- “@nonlfun”表示非线性部分的约束,我们也编写一个“m文件”
function [c,ceq] = nonlfun1(x)% 注意:这里的c实际上就是非线性不等式约束,ceq实际上就是非线性等式约束% 输入值x实际上就是决策变量,由x1和x2组成的一个向量% 返回值有两个,一个是非线性不等式约束c,一个是非线性等式约束ceq% nonlfun1是函数名称,到时候会被fmincon函数调用, 可以任意取名,但不能和目标函数fun1重名% 保存的m文件和函数名称得一致,也要为nonlfun1.mc = ceq = [];
end
①c,ceq中可能有多个约束,我们写成列向量的形式
②[x,fval] = fmincon(@fun1,x0,A,b,[],[],[][],@nonlfun1,option)
- 注意要把下标写在括号里
- 若不存在某种约束,则可用“[]”代替
4.基础例题
function [c,ceq] = nonlfun2(x)% 非线性不等式约束c = [-x(1)^2+x(2)-x(3)^2; % 一定要注意写法的规范,再次强调这里的x是一个向量!不能把x(1)写成x1x(1)+x(2)^2+x(3)^2-20];% 非线性等式约束ceq = [-x(1)-x(2)^2+2;x(2)+2*x(3)^2-3];
end
【1】默认方法求例1
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 -x2ceq = []; % 不存在非线性等式约束,所以用[]表示
end
function f = fun1(x)% 注意:这里的f实际上就是目标函数,函数的返回值也是f% 输入值x实际上就是决策变量,由x1和x2组成的向量% fun1是函数名称,到时候会被fmincon函数调用, 可以任意取名% 保存的m文件和函数名称得一致,也要为fun1.m
% max f(x) = x1^2 +x2^2 -x1*x2 -2x1 -5x2f = -x(1)^2-x(2)^2 +x(1)*x(2)+2*x(1)+5*x(2) ;
end
%% 例题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文件和nonlfun1.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)
% 结果也是可以计算出来的,但并不推荐这样做~
【2】用其他方法求例2
%% 使用其他算法对例题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.
% 这说明这个算法不适用我们这个约束条件,所以以后遇到了不能求解的情况,记得更换其他算法试试!!!
【3】改变初始值的影响
%% 生成不同的随机初始值来优化代码,有一定几率会触发上面那个Bug,因此不推荐
n = 10; % 重复n次
Fval = +inf; X = [0,0]; %初始化最优的结果
A = [-2 3]; b = 6;
for i = 1:nx0 = [rand()*10 , rand()*10]; %用随机数生成一个初始值(随机数的范围自己根据题目条件设置) [x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option); % 注意 fun1.m文件和nonlfun1.m文件都必须在当前文件夹目录下if fval < Fval % 如果找到了更小的值,那么就代替最优的结果Fval = fval;X = x;end
end
Fval = -Fval
X
【4】使用蒙特卡罗的方法来找初始值(推荐)
%% 使用蒙特卡罗的方法来找初始值(推荐)
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:nx = [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更新为初始值endend
end
disp('蒙特卡罗选取的初始值为:'); disp(x0)
A = [-2 3]; b = 6;
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1)
fval = -fval
【5】其他例题的求解
function f = fun2(x)% f = x(1)^2+x(2)^2 +x(3)^2+8 ; f = sum(x.*x) + 8; % 可别忘了x实际上是一个向量,我们可以使用矩阵的运算符号对其计算
end
function [c,ceq] = nonlfun2(x)% 非线性不等式约束c = [-x(1)^2+x(2)-x(3)^2; % 一定要注意写法的规范,再次强调这里的x是一个向量!不能把x(1)写成x1x(1)+x(2)^2+x(3)^2-20];% 非线性等式约束ceq = [-x(1)-x(2)^2+2;x(2)+2*x(3)^2-3];
end
%% 例题二的求解
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:nx = [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更新为初始值endend
end
disp('蒙特卡罗选取的初始值为:'); disp(x0)
lb = [0 0 0]; % 决策变量的下界
[x,fval] = fmincon(@fun2,x0,[],[],[],[],lb,[],@nonlfun2) % 注意 fun2.m文件和nonfun2.m文件都必须在当前文件夹目录下
function f = fun3(x)f = -prod(x); % 可别忘了x实际上是一个向量(prod表示连乘符号,用法和sum类似)
end
%% 例题三的求解(蒙特卡罗模拟那一讲的例题)
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
二.例题
1.选址问题
function f = fun5(xx) % 注意为了避免和下面的x同号,我们把决策变量的向量符号用xx表示(注意xx的长度为16)a=[1.25 8.75 0.5 5.75 3 7.25]; % 工地的横坐标b=[1.25 0.75 4.75 5 6.5 7.25]; % 工地的纵坐标x = [xx(13) xx(15)]; % 新料场的横坐标y = [xx(14) xx(16)]; % 新料场的纵坐标c = []; % 初始化用来保存工地和料场距离的向量 (这个向量就是我们的系数向量)for j =1:2for i = 1:6c = [c; sqrt( (a(i)-x(j))^2 + (b(i)-y(j))^2)]; % 每循环一次就在c的末尾插入新的元素endend% 下面我们要求吨千米数,注意c是列向量,我们计算非线性规划时给定的初始值x0是行向量f = xx(1:12) * c;
end
%% 选址问题
clear;clc
format long g %可以将Matlab的计算结果显示为一般的长数字格式(默认会保留四位小数,或使用科学计数法)
% % (1) 系数向量(原来线性规划问题的写法,我们只需要在此基础上改动一点就可以了)
% a=[1.25 8.75 0.5 5.75 3 7.25]; % 工地的横坐标
% b=[1.25 0.75 4.75 5 6.5 7.25]; % 工地的纵坐标
% x = [5 2]; % 料场的横坐标
% y = [1 7]; % 料场的纵坐标
% c = []; % 初始化用来保存工地和料场距离的向量 (这个向量就是我们的系数向量)
% for j =1:2
% for i = 1:6
% c = [c; sqrt( (a(i)-x(j))^2 + (b(i)-y(j))^2)]; % 每循环一次就在c的末尾插入新的元素
% end
% end
% (2) 不等式约束
A =zeros(2,16); % 注意这里要改成16
A(1,1:6) = 1;
A(2,7:12) = 1;
b = [20,20]';
% (3) 等式约束
Aeq = zeros(6,16); % 注意这里要改成16
for i = 1:6Aeq(i,i) = 1; Aeq(i,i+6) = 1;
end
beq = [3 5 4 7 6 11]'; % 每个工地的日需求量
%(4)上下界
lb = zeros(16,1);
% lb = [zeros(12,1); -inf*ones(4,1)]; 两个新料场坐标的下界可以设为-inf% 进行求解
% 注意哦,这里我们只尝试了这一个初始值,大家可以试试其他的初始值,有可能能够找到更好的解。
% 未来我会在遗传算法中再来看这个例题。
x0 = [3 5 0 7 0 1 0 0 4 0 6 10 5 1 2 7]; % 用第一问的结果作为初始值
[x,fval] = fmincon(@fun5,x0,A,b,Aeq,beq,lb) % 注意没有非线性约束,所以这里可以用[]替代,或者干脆不写
reshape(x(1:12),6,2) % 将x的前12个元素变为6行2列便于观察(reshape函数是按照列的顺序进行转换的,也就是第一列读完,读第二列,即x1对应x_1,1,x2对应x_2,1)
% 新坐标(5.74,4.99) (7.25,7.25)
% fval =
% 89.9231692432933
% 第一问的fval =
% 135.281541790676
135.281541790676 - 89.9231692432933 % 45.3583725473827
2.飞行管理问题
function [c,ceq] = nonlfun6(delta) % 决策变量delta为六架飞机调整的角度x = [150 85 150 145 130 0]; % 飞机初始位置的横坐标y = [140 85 155 50 150 0]; % 飞机初始位置的纵坐标theta = [243 236 220.5 159 230 52] * pi / 180; % 飞机初始的飞行方向角 v = 800; % 飞机速度co = cos(theta + delta); % 包含6个元素的向量si = sin(theta + delta); % 包含6个元素的向量% 下面开始计算飞机i和j之间的最短距离(只需要计算矩阵的一半即可)d = zeros(6); % 初始化飞机两两之间的最短距离矩阵for i = 2: 6for j = 1: i-1% 套用我们推导出来的公式计算飞机i和飞机j相距最近的时间fenzi = ((y(j)-y(i))*(si(j)-si(i)) +(x(j)-x(i))*(co(j)-co(i))) ; % 分子fenmu = v * ((co(j)-co(i))^2 + (si(j)-si(i))^2); % 分母t(i,j) =- fenzi / fenmu;if t(i,j) <0 d(i, j) = 1000; % 此时最初的位置就是相距最近的点,因为最初的时候所有飞机两两之间的距离就大于8,因此未来绝不会相撞,我们令它们的距离为一个特别大的数elsed(i, j) = sqrt((x(j)-x(i)+v*t(i,j)*(co(j)-co(i)))^2+(y(j)-y(i)+v*t(i,j)*(si(j)-si(i)))^2); end endend% 非线性不等式约束c =ones(15,1)*8.000001 - [d(2,1); d(3,1:2)'; d(4,1:3)'; d(5,1:4)'; d(6,1:5)']; % 12个非线性不等式约束: “最短距离>8” 等价于 “8 - 最短距离<0”% 注意: 由于Matlab标准型中取的是小于等于号,因此这里取一个比8略大的数:8.000001-最短距离<=0 ceq = []; % 没有非线性等式约束
end
function f = fun6(delta) % 决策变量delta为六架飞机调整的角度% f =sum(abs(delta)) * 180 / pi; % 目标函数第一种定义:绝对值的和(将弧度转换为度数)f = sum(delta .* delta) * (180 / pi)^2; % 目标函数第二种定义:平方和(将弧度转换为度数)
end
%% 飞行管理问题
format long g
%% (1)画六架飞机的位置
clear;clc
figure(1) % 生成一个图形
box on % 显示为封闭的盒子
% 绘制飞机的初始位置
data = [150 140 243;85 85 236;150 155 220.5;145 50 159;130 150 230;0 0 52];
plot(data(:,1),data(:,2),'.r')
axis([0 160,0,160]);% 设置坐标轴刻度范围
hold on;
% 在图上标上注释
for i = 1:6txt = ['飞机',num2str(i)];text(data(i,1)+2,data(i,2)+2,txt,'FontSize',8)
end
% 把Matlab做出来的图可以导出,然后再放到PPT中画出飞机飞行方向的箭头(就和讲义上的类似)%% 求解非线性规划问题
x0 = [0 0 0 0 0 0]; % 初始值
lb = -pi/6*ones(6,1);
ub = pi/6*ones(6,1);
[x,fval] = fmincon(@fun6,x0,[],[],[],[],lb,ub,@nonlfun6)
x = x * 180 / pi % 将弧度转换为度数
% 定义一:fval = 3.7315°
% 定义二: fval = 6.9547((°)^2)
数学建模更新12(非线性规划)相关推荐
- Python小白的数学建模课-12.非线性规划
非线性规划是指目标函数或约束条件中包含非线性函数的规划问题,实际就是非线性最优化问题. 从线性规划到非线性规划,不仅是数学方法的差异,更是解决问题的思想方法的转变. 非线性规划问题没有统一的通用方法, ...
- 建模计算机模拟框图,数学建模第12讲 计算机模拟.ppt
文档介绍: 数学建模与数学实验 后勤工程学院数学教研室 计算机模拟 些该胜卖咀合姓刽软脾姻闯晋乡昏缘易囊海屹盲救日斤状果揭减注变轰秦数学建模第12讲计算机模拟数学建模第12讲计算机模拟 实验目的 实验 ...
- 数学建模(12)——多元回归分析法
数学建模(12)--多元回归分析法 1.多元回归分析法的含义: 多元回归分析其实就是,抓重点.比如有各种各样的变量x1.x2.x3.x4,,,,要看看跟y有什么关系,多元回归就是把跟y相关的留下,无关 ...
- 数学建模常用算法—非线性规划
目录 模型的含义 问题转化为模型 模型建立 模型例题与求解 上篇我们学习了线性规划模型,这篇我们一起来学习一下非线性规划.非线性规划在工程.管理.经济.科研.军事等方面都有广泛的应用,为最优设计提供了 ...
- 【数学建模】8 非线性规划及例题讲解
1 定义 目标函数和约束条件中至少有一个非线性规划函数的数学规划问题称为非线性规划问题(UP问题) 2 MATLAB软件求解 函数 fmincon (1)建立M文件fun.m function f = ...
- 【数学建模】12 线性规划模型的求解方法
目录 1 图解法 2 MATLAB函数求解方法 3 Lingo法 4 课后习题 1 图解法 (1)有线性规划模型 • 目标函数 • 约束条件 在二元的约束条件画出来是直线,三元的约束条件画出来是一个平 ...
- 数学建模更新13(MATLAB绘制三维图【上】)
MATLAB绘制三维图 一.mesh函数以及拓展函数 1.mesh(X,Y,Z)的用法 [1]X是n维向量,Y是m维向量,Z是m*n维的矩阵 [2]X.Y和Z都是m*n维的矩阵 2.mesh(Z)的用 ...
- 数学建模之:非线性规划Python代码
from scipy import optimize as opt import numpy as np from scipy.optimize import minimize# 目标函数 def o ...
- 数学建模——更新1——excel直方图
目录 1.插入直方图 2.绘制频率分布直方图 步骤 1. 2. 3.生成表 4.生成频率 5.生成区间字符串 6.画柱状图 1.插入直方图 2.绘制频率分布直方图 步骤 1. 2. 3.生成表 4.生 ...
最新文章
- 需要反射时使用dynamic
- HTML5 监听当前位置
- 山西计算机软考知识点,计算机软考考试必备知识点:数据标准化
- 远程桌面中Tab键不能补全的解决办法
- 深度学习激活函数比较
- 零基础可以学python吗-学Python需要什么基础知识?零基础可以学Python吗?
- attachment old API read - DB debug
- php mysql 绕过_PHP中md5绕过
- 吊打6599元的三星?买手机莫慌 三款国产新手机将发
- CEPH RGW集群和bucket的zone group 不一致导致的404异常解决 及 使用radosgw-admin metadata 命令设置bucket metadata 的方法
- 了解不同种类的windows存储驱动
- CAM350 12.1免费下载
- 将输入的字符串逆序输出
- 用JavaScript写一个简单的网页倒计时插件
- Mybatis Plus分页Page total始终为0
- IBM AIX初级培训总结
- 苹果支付IAP V1
- 锂电池放空后充不进电_锂电池放置太久无法充电 血的教训!
- 计算机硬件的维护知识,计算机硬件维护的知识
- 一台台式计算机应该具有哪些设备,电脑硬件有哪些?组装一台电脑需要哪些配件详解...