数模3—Matlab线性规划、非线性规划、多目标规划（超全解法合集）

线性规划，非线性规划，多目标规划都归于优化类模型

一、线性规划

🎐例题

张麻子既要攻碉楼又要追替身，他们一伙6人，总共1200发子弹;每有一人攻碉楼会给百姓带来40点士气值，每有一人追替身会给百姓带来30点士气值;攻碉楼每人需240发子弹，追替身每人需120发。
问攻碉楼和追替身各派几个人，能使百姓的士气值最大?

解析为：

注意：
线性规划就是在一组线性约束条件下,求线性目标函数的最大或最小值
"线性"意味着所有变量都是一次方

🎐线性规划适用的典型赛题

题目中提到"怎样安排/分配”"尽量多(少)” “最多(少)” "利润最大” “最合理” 等词;

生产安排：原材料、设备有限制，总利润最大
若生产两种机床，利润分别为XXX; A机器和B机器加工，有顺序要求，有不同损耗费用，不同的工作时间…;问题：怎样安排生产使得总利润最大?

投资收益:资产配置、组合投资、总收益最大（收益率、损失率可能不是线性了）
若总资金为M，有n种资产可以配置。每种资产的平均收益率…，风险损失率…，手续费…;
问题:设计组合投资方案，使得收益尽可能大(本质是多目标规划，可化简为一个目标的线性规划)

销售运输:产地、销地、产量、销量、运费，总运费最省
商品有m个产地和n个销地，需要从产地运到销地。各产地的产量…，各销地需求量…，由a产地运到b销地的运价xxx;问题:如何调运才能使总运费最省?

车辆安排:路线、起点终点、承载量、时间点、车次安排最合理（合理就是指最多/最少！）
不同种类的车辆有各自的承载量，工地各点之间要安排车辆运输。工地里有多条路线…满足用工需求的情况下…;问题:如何安排车辆能使产量尽可能大?

(注意:一般资料里所说的整数规划、0-1规划往往也默认为线性规划的特例)

🎐如何用Matlab求解线性规划

模型MATLAB求解: Linprog函数 注意：只能求最小值！
PS：matlab求线性规划要先把模型，弄出标准型。所以！线性规划嘛，都用lingo解。

求最大值：y的最大值等价于-y的最小值
x>=b,也等价于-x<=b
[x, fval] = linprog (f,A, b, Aeq, beq, Ib, ub)
f：目标函数的系数列向量
f = [-40; -30]
A:约束条件里的变量系数矩阵
b:约束条件里的常数量矩阵
Aeq：等式约束条件
beq：等式约束条件
lb：
ub： 已经在不等式约束里面写了就可以不写。

最后的到的x是一个矩阵。

---

✨若不存在不等式约束，用"[]”代替A和b: [x, fval ] = linprog (f,[],[], Aeq, beq, lb, ub)
✨若不存在等式约束，用"[]”代替Aeq和beq: [x, fval ] = linprog Cf,A,b,[],[] ,lb,ub)
✨没有等式约束和最小、最大取值的约束时，可以不写Aeq, beq和lb, ub: [x, fval ] = linprog (f,A,b)

✨若题目求最大值:目标函数等号两端加符号转为求最小值;求解后目标值再取负！！千万别忘记

代码：

f = [-40;-30]; %目标函数中变量的系数矩阵
a = [1,1;-1,0;0,-1;240, 120];
%小于等于的约束条
b = [6;-1;-1;1200];
%小于等于的约束条
%标准形式: [x，fval ]= linprog (f,A,b, aeq, beq,1b
%从某个变量开始后面都没有?全都不写
%没有等式约束条件、但有上下界约束时，用[]在函数里占
[x,y] = linprog(f,a,b)
y=-y
% zeros(m, n)是元素全为0的m行n列的矩阵
1
2
3
4
5
6
7
8
9
10
11

结果：

x =
    4.0000
    2.0000
y =
  -220
y =
   220
1
2
3
4
5
6
7

结论，4个人攻碉楼，2个人追替身的情况下，百姓的士气最大，最大值为220。

🎐lingo求解线性规划

直接把建模的式子，搬过来。注意！乘号一定要打出来。
分号也一定要一个式子一个

max = 40*x1 + 30*x2;
x1 + x2 <= 6;
x1 >= 1;
x2 >=1;
240*x1 + 120*x2 <=1200;
1
2
3
4
5

结果如图：

二、非线性规划

和线性规划一样都是求最大的收益

线性规划:所有变量都是一次方
非线性规划:至少一个变量不是一次方

🎐线性规划适用的典型赛题

题目中提到"怎样安排/分配”"尽量多(少)” “最多(少)” "利润最大” “最合理” 等词;但变量非一次方

生产安排:原材料、设备有限制，总利润最大(目标函数或约束条件非线性)
生产两种机床，利润分别为XXX, A机器和B机器加工，两种机器工作时间…;成本或利润与某变量的关系是非线性的，例如A机器每小时生产费用是f= 50x + 0.2x2 (元)，求安排生产使得总利润最大

投资规划:资产配置、收益率、损失率、组合投资、总收益率最大\最佳投资方案
总资金为M，有n种资产可以配置，平均收益率…，风险损失率…，手续费…，设计组合投资方案使得收益尽可能大，总体风险尽可能小

角度调整:飞行管理避免相撞; 影院最佳视角
飞机位置，速度，进入区域后判定是否相撞， 飞机飞行方向角调整的幅度尽量小
电影院视角、仰角影响观影体验，什么位置观影最佳涉及三角函数，为非线性

🎐lingo求解非线性规划

lingo代码：

min = x1^2 + x2^2 + x3^2 +8;
x1^2 - x2 + x3^2 >=0;
x1 + x2^2 + x3^3 <=20;
-x1 - x2^2 + 2=0;
x2 + 2*x3^2 =3;
x1>=0;x2>=0;x3>=0;
1
2
3
4
5
6

结果如下：

🎐Matlab求解

注意：即是初始值不一样，反正都是遍历迭代查找，所以最终结果都是一样的。
[x, fval] = fmincon( fun,x0, A, b, Aeq, beq, lb, ub, nonlcon)

• fun 单独脚本文件里定义的目标函数
• x0 决策变量的初始值(不影响最终结果)
• A, b 线性约束的不等式变量系数矩阵和常数项矩阵(<=或<)
• Aeq, beq 线性约束的等式变量系数矩阵和常数项矩阵
• lb, ub 决策变量的最小取值和最大取值
• nonlcon 非线性约束，包括不等式和等式

fun1 目标函数

function f = fun1(x)
f = sum(x.^2) + 8;
%目标函数
1
2
3

非线性约束条件

function [g,h]=fun2(x)
g = [-x(1)^2+x(2)-x(3)^2
x(1)+x(2)^2+x(3)^3-20]    %g非线性不等式

h=[-x(1)-x(2)^2+2
x(2)+2*x(3)^2-3]       %h是非线性等式
1
2
3
4
5
6

这个例子没有线性约束条件

求解代码：

[x,y] = fmincon('fun1',[0;0;0],[],[],[],[],[0;0;0],[],'fun2')
1

结果：

x =
    0.5522
    1.2033
    0.9478
y =
   10.6511
1
2
3
4
5
6

注意： 那些系数，都是列变量哦！并且：

lingo和Matlab次方都是一样的表示。但是在变量表示上Matlab略胜一筹啦~~因为Matlab的向量的表示十分丝滑，矩阵运算十分完善。

三、多目标规划

本质:既要XXX，又要XXX

🎐多目标规划基础实例

模型解析：
多目标:现在设有3个目标:
●1.尽量使产品I的产量不超过产品II的产量;
●2.尽可能充分利用设备，但不希望加班
●3.尽可能使利润不少于56万

上面这张图的第三个写错了，应该是P3d3-
查文献或者自己主观设优先因子，说是查了文献的

🎐多目标规划适用的典型赛题

这个类型建模比较重要。求解可以用Matlab，也可以用lingo。求解办法和上面的线性规划和非线性规划类似。https://shenlong.blog.csdn.net/article/details/110501496
建议使用数学规划法，求出具体答案！

也可以！！ 直接用Matlab的遗传算法求解：参考下面这篇文章
https://zstar.blog.csdn.net/article/details/119302954

🎐基于NSGA-II算法的gamultiobj函数求解多目标优化

1.多目标规划的问题描述：

2.求解方法
Matlab中提供函数gamultiobj采用的算法就是基于NSGA-II改进的一种多目标优化算法（a variant of NSGA-II）,接下来对目标规划中的一些概念进行介绍。

目前求解多目标优化问题方法算法主要有基于数学的规划方法和基于遗传算法的两类方法；其中带精英策略的快速非支配排序算法（NSGA-II)是影响最大和应用范围最广的一种多目标遗传算法。算法主要优点：

1.提出了快速非支配的排序算法，降低了计算非支配序的复杂度。
2.引入了精英策略，扩大了采样空间。将父代种群与其产生的子代种群组合在一起，共同通过竞争来产生下一代种群，这有利于是父代中的优良个体得以保持，保证那些优良的个体在进化过程中不被丢弃，从而提高优化结果的准确度。并且通过对种群所有个体分层存放，使得最佳个体不会丢失，能够迅速提高种群水平。
3.引入拥挤度和拥挤度比较算子，这不但克服了NSGA算法中需要人为指定共享参数的缺陷，而且将拥挤度作为种群中个体之间的比较准则，使得准Pareto域中的种群个体能均匀扩展到整个Pareto域，从而保证了种群的多样性。

3.求解案例
求解的通用形式式

例1
使用gamultiobj求解,主要分为两部分：
Fun.m文件

%%定义求解函数文件
function y=Fun(x) 
y(1)=x(1)^4-10*x(1)^2+x(1)*x(2)+x(2)^4-x(1)^2*x(2)^2; 
y(2)=x(2)^4-x(1)^2*x(2)^2+x(1)^4+x(1)*x(2); 

1
2
3
4
5

main.m文件

%%求解文件
clear;clc;

fitnessfcn=@Fun;%适应度函数句柄
 nvars=2; %变量个数
 lb=[-5,-5]; %下限
 ub=[5,5]; %上限
 A=[];b=[]; %线性不等式约束
 Aeq=[];beq=[]; %线性等式约束
 
options=gaoptimset('paretoFraction',0.3,'populationsize',200,'generations',200,'stallGenLimit',200,'TolFun',1e-10,'PlotFcns',@gaplotpareto); 
% 最优个体系数paretoFraction为0.3;
%种群大小populationsize为100，
%最大进化代数generations为200，
% 停止代数stallGenLimit为200，
%适应度函数偏差TolFun设为1e-10，
%函数gaplotpareto：绘制Pareto前沿
[x,fval]=gamultiobj(fitnessfcn,nvars,A,b,Aeq,beq,lb,ub,options)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

句柄的概念：https://blog.csdn.net/xiaofeiyu321/article/details/89791451
求解结果：

图像会实时记录解的情况。

最后的求解结果，即帕累托最优解集储存在[x,fval]中，fval是x对应的目标函数值。fval大致构成了一条空间曲线——Pareto前沿。若各个解较为均匀分布，说明该图包含了大部分最优解情况，全局性优，适用性强。
在满足Pareto最优的条件下，是没有办法在不让某一优化目标受损的情况下，令另一方目标获得更优的。所以这些解均为最优，对最优解的具体选择可以根据实际情况。

更多求解案例，可参考：

https://zstar.blog.csdn.net/article/details/119302954

如果是求最大，可参考这个的第一个例子;
https://blog.csdn.net/weixin_47102975/article/details/108253584

这两个和主函数放在一个文件，在主函数的前面

%% 模型设置
% 适应度函数的函数句柄
fitnessfcn=@Fun;
% 变量个数
nvars=4;
% 约束条件形式1：下限与上限（若无取空数组[]）
% lb<= X <= ub
lb=[0,0,0,0];
ub=[];

% 约束条件形式2：线性不等式约束（若无取空数组[]）
% A*X <= b 
A = [0    0 1 1
    -1/3  0 0 0
     0 -1/2 0 0
     0    0 0 0];

b = [48 ; 30 ; 30 ; 0];

% 约束条件形式3：线性等式约束（若无取空数组[]）
% Aeq*X == beq
Aeq=[1 1 0 0;0 0 0 0; 0 0 0 0; 0 0 0 0];
beq=[120;0;0;0];
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

求解器设置部分：

%% 求解器设置
% 最优个体系数paretoFraction
% 种群大小populationsize
% 最大进化代数generations
% 停止代数stallGenLimit
% 适应度函数偏差TolFun
% 函数gaplotpareto：绘制Pareto前沿 
options=gaoptimset('paretoFraction',0.3,'populationsize',200,'generations',300,'stallGenLimit',200,'TolFun',1e-10,'PlotFcns',@gaplotpareto);
1
2
3
4
5
6
7
8

gamultiobj求解与结果输出部分：逆向重点！！

%% 主求解
[x,fval]=gamultiobj(fitnessfcn,nvars,A,b,Aeq,beq,lb,ub,options)
%% 结果提取
% 因为gamultiobj是以目标函数分量取极小值为目标，
% 因此在y=Fun(x)里取相反数的目标函数再取相反数画出原始情况
plot(-fval(:,1),fval(:,2),'pr')
xlabel('f_1(x)')
ylabel('f_2(x)')
title('Pareto front')
grid on
1
2
3
4
5
6
7
8
9
10

目标函数: (这一段需放在脚本最后或单独放在一个文件里)

function y=Fun(x)
	% y是目标函数向量。有几个目标函数y就有多少个维度（数组y的长度）
	% 因为gamultiobj是以目标函数分量取极小值为目标，
	% 因此有些取极大值的目标函数注意取相反数
	y(1)=-(x(1)*100/3 + x(3)*90/3  + x(2)*80/2+x(4)*70/2);
	y(2)=x(3)+x(4);
end
1
2
3
4
5
6
7