数学规划模型

￡神魔★判官ぃ 2023-09-28 16:28 18阅读 0赞 
  1. 来源：数学建模清风学习内容整理

文章目录

  • 一、概述
      • （1）什么是数学规划？
      • （2）数学规划的一般形式
      • （3）数学规划的分类
  • 二、线性规划问题的求解
      • （1）Matlab中线性规划的标准型
      • （2）Matlab求解线性规划的函数
    • 三、线性规划的典型例题
      • （1）例题一：生产决策问题
      • （2）例题二：投料问题
  • 四、非线性规划问题的求解
      • （1）Matlab中非线性规划的标准型
      • （2）Matlab中求解非线性规划的命令
            • 对于例题的matlab使用讲解
    • 五、非线性规划的典型例题
      • （1）例题一：选址问题
      • （2）例题二：飞行管理问题
  • 六、整数规划
      • （1）Matlab中线性整数规划求解
      • （2）Matlab中线性 0 - 1 规划求解
    • 七、整数规划的典型例题
      • （1）例题一：背包问题
      • （2）例题二：指派问题
      • （3）例题三：钢管切割问题
  • 八、最大化最小化模型
      • （1）模型的一般形式
      • （2）典型例题
      • （3）模型的求解
  • 九、多目标规划模型
      • （1）求解思路
      • （2）典型例题
  • 十、拓展例题
      • （1）护士值班问题
      • （2）货舱装运问题
      • （3）非线性规划的求解
      • （4）覆盖问题

一、概述

（1）什么是数学规划？

  • 求目标函数在一定约束条件下的极值问题

数学规划是运筹学的一个分支，用来研究：在给定的条件下（约束条件），如何按照某一衡量指标（目标函数）来寻求计划、管理工作中的最优方案。

（2）数学规划的一般形式

在这里插入图片描述

（3）数学规划的分类

① 线性规划（Linear programming）
如果目标函数 f ( x ) f(x) f(x)和约束条件均是决策变量的线性表达式，那么此时的数学规划问题就属于线性规划。

② 非线性规划（nonlinear programming）
当目标函数 f ( x ) f(x) f(x)或者约束条件中有一个是决策变量 x x x的非线性表达式，那么此时的数学规划问题就属于非线性规划。

③ 整数规划（integer programming）
整数规划是一类要求变量取整数值的数学规划（线性整数规划和非线性整数规划）。

④ 0-1规划（0-1 programming）
整数规划的特例，整数变量的取值只能为0和1。

二、线性规划问题的求解

（1）Matlab中线性规划的标准型

在这里插入图片描述
在这里插入图片描述

  1. % 1题代码
  2. [x,fval]=linprog(c,A,b,[],[],lb)    %不写ub意味着不存在上界约束
  4. % 2题代码
  5. [x,fval]=linprog(c,A,b,Aeq,beg,lb)
  7. % 3题代码
  8. [x,fval]=linprog(c,A,b,Aeq,beg,lb)
  9. fval=-fval

（2）Matlab求解线性规划的函数

  1. [x,fval]=linprog(c,A,b,Aeq,beq,lb,ub,xo)

① xo表示给定Matlab迭代求解的初始值（一般不用给）；
② c，A，b，Aeq，beq，lb，ub的意义和标准型中的意义一致 ；
③ 若不存在不等式约束，可用 [ ] 替代 A 和 b ；
④ 若不存在等式约束，可用 [ ] 替代 Aeq 和 beq ；
⑤ 若某个 x i x_i xi无下届或上届，则设置 lb(i)=-inf，ub(i)=+inf ；
⑥ 返回的 x x x表示最小值处的 x x x取值；fval 表示最优解处时取得的最小值 ；
⑦ 不是所有的线性规划问题都有唯一解，可能无解或有无穷多的解 ；
⑧ 如果求的是最大值，别忘了在最后给 fval 加一个负号 ；

三、线性规划的典型例题

（1）例题一：生产决策问题

在这里插入图片描述

  1. %% 生产决策问题
  2. format long g   %可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）
  3. % (1) 系数向量
  4. c = zeros(9,1); % 初始化目标函数的系数向量全为0
  5. c(1) = 1.25 -0.25 -300/6000*5;  % x1前面的系数是c1
  6. c(2) = 1.25 -0.25 -321/10000*7;
  7. c(3) = -250 / 4000 * 6;
  8. c(4)  = -783/7000*4;
  9. c(5) = -200/4000 * 7;
  10. c(6) = -300/6000*10;
  11. c(7) = -321 / 10000 * 9;
  12. c(8) = 2-0.35-250/4000*8;
  13. c(9) = 2.8-0.5-321/10000*12-783/7000*11;
  14. c = -c;  % 我们求的是最大值，所以这里需要改变符号
  15. % (2) 不等式约束
  16. A = zeros(5,9);
  17. A(1,1) = 5;  A(1,6) = 10;
  18. A(2,2) = 7;  A(2,7) = 9; A(2,9) = 12;
  19. A(3,3) = 6;  A(3,8) = 8;
  20. A(4,4) = 4;  A(4,9) = 11;
  21. A(5,5) = 7;  
  22. b = [6000 10000 4000 7000 4000]';
  23. % (3) 等式约束
  24. Aeq = [1 1 -1 -1 -1 0 0 0 0;
  25.             0 0 0 0 0 1 1 -1 0];
  26. beq = [0 0]';
  27. %（4）上下界
  28. lb = zeros(9,1);
  30. % 进行求解
  31. [x fval] = linprog(c, A, b, Aeq, beq, lb)
  32. fval = -fval
  33. % fval =
  34. %           1146.56650246305
  35. %  注意，本题应该是一个整数规划的例子，我们在后面的整数规划部分再来重新求解。
  36. intcon = 1:9;
  37. [x,fval]=intlinprog(c,intcon,A,b,Aeq,beq,lb)
  38. fval = -fval

（2）例题二：投料问题

在这里插入图片描述

在这里插入图片描述

  1. %% 投料问题
  2. clear,clc
  3. format long g   %可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）
  4. % (1) 系数向量
  5. a=[1.25  8.75  0.5  5.75  3  7.25];  % 工地的横坐标
  6. b=[1.25  0.75  4.75    5  6.5  7.25];   % 工地的纵坐标
  7. x = [5  2];  % 料场的横坐标
  8. y = [1  7];  % 料场的纵坐标
  9. c = [];  % 初始化用来保存工地和料场距离的向量 (这个向量就是我们的系数向量）
  10. for  j =1:2
  11.     for i = 1:6
  12.         c = [c;  sqrt( (a(i)-x(j))^2 + (b(i)-y(j))^2)];  % 每循环一次就在c的末尾插入新的元素
  13.     end
  14. end
  15. % (2) 不等式约束
  16. A =zeros(2,12);
  17. A(1,1:6) = 1;
  18. A(2,7:12) = 1;
  19. b = [20,20]';
  20. % (3) 等式约束
  21. Aeq = zeros(6,12);  
  22. for i = 1:6
  23.     Aeq(i,i) = 1;  Aeq(i,i+6) = 1;
  24. end
  25. % Aeq = [eye(6),eye(6)]  % 两个单位矩阵横着拼起来
  26. beq = [3 5 4 7 6 11]';  % 每个工地的日需求量
  27. %（4）上下界
  28. lb = zeros(12,1);
  30. % 进行求解
  31. [x fval] = linprog(c, A, b, Aeq, beq, lb)
  32. x = reshape(x,6,2)  % 将x变为6行2列便于观察（reshape函数是按照列的顺序进行转换的，也就是第一列读完，读第二列，即x1对应x_1,1，x2对应x_2,1）
  34. % fval =
  35. %           135.281541790676

四、非线性规划问题的求解

（1）Matlab中非线性规划的标准型

在这里插入图片描述

在这里插入图片描述

（2）Matlab中求解非线性规划的命令

  1. [x,fval]=fmincon(@fun,x0,A,b,Aeq,beq,lb,ub,@nonlfun,option)

（1） 非线性规划中对于初始值 x 0 x_0 x0 的选取非常重要，因为非线性规划的算法求解出来的是一个局部最优解。
（线性规划不存在这个问题）
（2） 如果要求“全局最优解”，有两种思路：
① 给定不同的初始值，在里面找到最优解；
② 先用蒙特卡罗模拟，得到一个蒙特卡罗解，再将这个解作为初始值来求最优解；
（3）“option” 选项可以给定求解的算法，一共有四种：
interior-point（内点法）
sqp（序列二次规划法）
active-set（有效集法）
trust-region-reflective（信赖域反射算法）
（4） 不同的算法有其各自的优缺点和适用情况，我们可以改变求解的算法来看求解的结果是否变好；
（5 ）“ @fun ” 表示目标函数，要编写一个独立的“.m”文件储存目标函数：

  1. funciton f=fun(x)       % 注：1.fun可以任意取名；
  2.          f= ...         %    2.f也可以任意取名，但返回的f和函数内部的f得完全一致；
  3. end                     %    3.这里的x实际上是表示决策变量的向量，其行列方向取决于初始值x；
  4.                         %    4.调用函数：fmincon(@fun, ... ) 求解

（6）“@nonlfun”表示非线性部分的约束，同样得编写一个独立的 “.m” 文件储存非线性约束条件：

  1. function[c,ceq]=nonlfun(x)             %  1. nonlfun同样可任意取名，不和上面的fun相同就可，保存的.m文件也得是这个名字；
  2.         c=[非线性不等式约束1；          %  2. c、ceq中可能有多个约束，因此写成列向量的形式；
  3.            非线性不等式约束p；]         %  3. 若不存在非线性不等式约束，则可令 c=[];
  4.       ceq=[非线性等式约束1；            %  4. 调用函数：fmincon（...,@nonlfun,options）求解
  5.            非线性等式约束q；         
  6. end

（7） 注意要把下标改写为括号，例如： f = x 1 2 + 3 x 2 f=x_1^2+3x_2 f=x12+3x2 写成matlab能识别的形式：f=x(1)^2+3*x(2)；
（8） 若不存在某种约束，则可用 “ [ ] ” 替代，若后面全为 “ [ ] ” 且不指定option（使用默认的求解方法），则 “ [ ] ” 也可以省略掉；

对于例题的matlab使用讲解

  • 主函数

    %% 非线性规划的函数
    % [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 -x1x2 -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)
    % 结果也是可以计算出来的，但并不推荐这样做~

    %% 使用其他算法对例题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

    1. x0 = [rand()*10 , rand()*10];  %用随机数生成一个初始值（随机数的范围自己根据题目条件设置） 
    2. [x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option); % 注意 fun1.m文件和nonlfun1.m文件都必须在当前文件夹目录下
    3. if fval < Fval  % 如果找到了更小的值，那么就代替最优的结果
    4.     Fval = fval;
    5.     X = x;
    6. 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

    1. x = [x1(i), x2(i)];  %构造x向量, 这里千万别写成了：x =[x1, x2]
    2. if ((x(1)-1)^2-x(2)<=0)  & (-2*x(1)+3*x(2)-6 <= 0)     % 判断是否满足条件
    3.     result = -x(1)^2-x(2)^2 +x(1)*x(2)+2*x(1)+5*x(2) ;  % 如果满足条件就计算函数值
    4.     if  result  < fmin  % 如果这个函数值小于我们之前计算出来的最小值
    5.         fmin = result;  % 那么就更新这个函数值为新的最小值
    6.         x0 = x;  % 并且将此时的x1 x2更新为初始值
    7.     end
    8. end

    end
    disp(‘蒙特卡罗选取的初始值为：’); disp(x0)
    A = [-2 3]; b = 6;
    [x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1)
    fval = -fval 

  1. %% 例题二的求解
  2. x0 = [1 1 1];  %任意给定一个初始值 
  3. lb = [0 0 0];  % 决策变量的下界
  4. [x,fval] = fmincon(@fun2,x0,[],[],[],[],lb,[],@nonlfun2)  % 注意 fun2.m文件和nonfun2.m文件都必须在当前文件夹目录下
  5. % x =
  6. %          0.552167405729277          1.20325915507969         0.947824046150443
  7. % fval =
  8. %           10.6510918606939
  12. %% 使用蒙特卡罗的方法来找初始值(推荐）
  13. clc,clear;
  14. n=1000000; %生成的随机数组数
  15. x1= unifrnd(0,2,n,1);   % 生成在[0,2]之间均匀分布的随机数组成的n行1列的向量构成x1
  16. x2 = sqrt(2-x1);  % 根据非线性等式约束用x1计算出x2
  17. x3 = sqrt((3-x2)/2); % 根据非线性等式约束用x2计算出x3
  18. fmin=+inf; % 初始化函数f的最小值为正无穷（后续只要找到一个比它小的我们就对其更新）
  19. for i=1:n
  20.     x = [x1(i), x2(i), x3(i)];  %构造x向量, 这里千万别写成了：x =[x1, x2, x3]
  21.     if (-x(1)^2+x(2)-x(3)^2<=0) & (x(1)+x(2)^2+x(3)^2-20<=0)   % 判断是否满足条件
  22.         result =sum(x.*x) + 8 ;  % 如果满足条件就计算函数值
  23.         if  result  < fmin  % 如果这个函数值小于我们之前计算出来的最小值
  24.             fmin = result;  % 那么就更新这个函数值为新的最小值
  25.             x0 = x;  % 并且将此时的x1 x2 x3更新为初始值
  26.         end
  27.     end
  28. end
  29. disp('蒙特卡罗选取的初始值为：'); disp(x0)
  30. lb = [0 0 0];  % 决策变量的下界
  31. [x,fval] = fmincon(@fun2,x0,[],[],[],[],lb,[],@nonlfun2)  % 注意 fun2.m文件和nonfun2.m文件都必须在当前文件夹目录下
  33. %% 例题三的求解(蒙特卡罗模拟那一讲的例题)
  34. clear;clc
  35. % 蒙特卡罗模拟得到的最大值为3445.6014
  36. % 最大值处x1 x2 x3的取值为：
  37. %           22.5823101903968          12.5823101903968          12.1265223966757
  38. A = [1 -2 -2;  1 2 2];  b = [0 72];
  39. x0 = [ 22.58   12.58  12.13];
  40. Aeq = [1 -1 0]; beq = 10;
  41. lb = [-inf 10 -inf];  ub = [inf 20 inf];  
  42. [x,fval] = fmincon(@fun3,x0,A,b,Aeq,beq,lb,ub,[])  % 注意没有非线性约束，所以这里可以用[]替代，或者干脆不写
  43. fval = -fval
  46. % % 注意：代码文件仅供参考，一定不要直接用于自己的数模论文中
  47. % % 国赛对于论文的查重要求非常严格，代码雷同也算作抄袭

  • 被调用的函数：fun1 （例题一的目标函数）

    function f = fun1(x)

    1. % 注意：这里的f实际上就是目标函数，函数的返回值也是f
    2. % 输入值x实际上就是决策变量，由x1和x2组成的向量
    3. % fun1是函数名称，到时候会被fmincon函数调用, 可以任意取名
    4. % 保存的m文件和函数名称得一致，也要为fun1.m

    % max f(x) = x1^2 +x2^2 -x1*x2 -2x1 -5x2

    1. f = -x(1)^2-x(2)^2 +x(1)*x(2)+2*x(1)+5*x(2) ;

    end

  • 被调用的函数：fun2 （例题二的目标函数）

    function f = fun2(x)

    1. %     f = x(1)^2+x(2)^2 +x(3)^2+8 ; 
    2. f = sum(x.*x) + 8;  % 可别忘了x实际上是一个向量，我们可以使用矩阵的运算符号对其计算

    end

  • 被调用的函数：fun3 （例题三的目标函数）

    function f = fun3(x)

    1. f = -prod(x);  % 可别忘了x实际上是一个向量（prod表示连乘符号，用法和sum类似）

    end

  • 被调用的函数： nonlfun1 （例题一的非线性约束条件）

    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 （例题二的非线性约束条件）

    function [c,ceq] = nonlfun2(x)

    1. % 非线性不等式约束
    2. c = [-x(1)^2+x(2)-x(3)^2;   % 一定要注意写法的规范，再次强调这里的x是一个向量！不能把x(1)写成x1
    3.         x(1)+x(2)^2+x(3)^2-20];
    4. % 非线性等式约束
    5. ceq = [-x(1)-x(2)^2+2;
    6.             x(2)+2*x(3)^2-3];

    end

五、非线性规划的典型例题

（1）例题一：选址问题

在这里插入图片描述
在这里插入图片描述

  • 主函数

    %% 选址问题
    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:6

    1. Aeq(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

  • 被调用的函数

    function f = fun5(xx) % 注意为了避免和下面的x同号，我们把决策变量的向量符号用xx表示（注意xx的长度为16）

    1. a=[1.25  8.75  0.5  5.75  3  7.25];  % 工地的横坐标
    2. b=[1.25  0.75  4.75    5  6.5  7.25];   % 工地的纵坐标
    3. x = [xx(13)  xx(15)];  % 新料场的横坐标
    4. y = [xx(14)  xx(16)];  % 新料场的纵坐标
    5. c = [];  % 初始化用来保存工地和料场距离的向量 (这个向量就是我们的系数向量）
    6. for  j =1:2
    7.     for i = 1:6
    8.         c = [c;  sqrt( (a(i)-x(j))^2 + (b(i)-y(j))^2)];  % 每循环一次就在c的末尾插入新的元素
    9.     end
    10. end
    11. % 下面我们要求吨千米数，注意c是列向量，我们计算非线性规划时给定的初始值x0是行向量
    12. f = xx(1:12) * c;

    end

（2）例题二：飞行管理问题

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述
计算新坐标的示意图

  • 主函数

    %% 飞行管理问题
    format long g
    %% （1）画六架飞机的位置
    clear;clc
    figure(1) % 生成一个图形
    box on % 显示为封闭的盒子
    % 绘制飞机的初始位置
    data = [150 140 243;

    1. 85    85    236;
    2. 150    155    220.5;
    3. 145    50    159;
    4. 130    150    230;
    5. 0    0    52];

    plot(data(:,1),data(:,2),’.r’)
    axis([0 160,0,160]);% 设置坐标轴刻度范围
    hold on;
    % 在图上标上注释
    for i = 1:6

    1. txt = ['飞机',num2str(i)];
    2. text(data(i,1)+2,data(i,2)+2,txt,'FontSize',8)

    end
    % 把Matlab做出来的图可以导出，然后再放到PPT中画出飞机飞行方向的箭头

    %% 求解非线性规划问题
    x0 = [0 0 0 0 0 0]; % 初始值
    lb = -pi/6ones(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)

  • 被调用函数：fun6

    function f = fun6(delta) % 决策变量delta为六架飞机调整的角度
    % f =sum(abs(delta)) 180 / pi; % 目标函数第一种定义：绝对值的和（将弧度转换为度数）
    f = sum(delta .     delta) * (180 / pi)^2; % 目标函数第二种定义：平方和（将弧度转换为度数）
    end

  • 被调用函数：nonlfun6

    function [c,ceq] = nonlfun6(delta) % 决策变量delta为六架飞机调整的角度

    1. x = [150 85 150 145 130  0]; % 飞机初始位置的横坐标
    2. y = [140 85 155  50 150  0]; % 飞机初始位置的纵坐标
    3. theta = [243 236 220.5 159 230 52] * pi / 180; % 飞机初始的飞行方向角 
    4. v = 800;  % 飞机速度
    5. co = cos(theta + delta);  % 包含6个元素的向量
    6. si = sin(theta + delta);  % 包含6个元素的向量
    7. % 下面开始计算飞机i和j之间的最短距离（只需要计算矩阵的一半即可）
    8. d = zeros(6);  % 初始化飞机两两之间的最短距离矩阵
    9. for i = 2: 6
    10.     for j = 1: i-1
    11.         % 套用我们推导出来的公式计算飞机i和飞机j相距最近的时间
    12.         fenzi = ((y(j)-y(i))*(si(j)-si(i)) +(x(j)-x(i))*(co(j)-co(i))) ;  % 分子
    13.         fenmu =  v * ((co(j)-co(i))^2 + (si(j)-si(i))^2);  % 分母
    14.         t(i,j) =- fenzi / fenmu;
    15.         if t(i,j) <0  
    16.             d(i, j) = 1000; % 此时最初的位置就是相距最近的点，因为最初的时候所有飞机两两之间的距离就大于8，因此未来绝不会相撞，我们令它们的距离为一个特别大的数
    17.         else
    18.             d(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); 
    19.         end 
    20.     end
    21. end
    22. % 非线性不等式约束
    23. 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)'];  
    24. % 12个非线性不等式约束: “最短距离>8” 等价于 “8 - 最短距离<0”
    25. % 注意： 由于Matlab标准型中取的是小于等于号，因此这里取一个比8略大的数：8.000001-最短距离<=0 
    26. ceq = [];  % 没有非线性等式约束

    end

六、整数规划

在这里插入图片描述

（1）Matlab中线性整数规划求解

  1. [x,fval]=intlinprog(c,intcon,A,b,Aeq,beq,lb,ub)

在这里插入图片描述

在这里插入图片描述

  1. %% 线性整数规划问题
  2. %% 例1
  3. c=[-20,-10]';
  4. intcon=[1,2];  % x1和x2限定为整数
  5. A=[5,4;
  6.       2,5];
  7. b=[24;13];
  8. lb=zeros(2,1);  
  9. [x,fval]=intlinprog(c,intcon,A,b,[],[],lb)
  10. fval = -fval
  12. %% 例2
  13. c=[18,23,5]';
  14. intcon=3;  % x3限定为整数
  15. A=[107,500,0;
  16.       72,121,65;
  17.       -107,-500,0;
  18.       -72,-121,-65];
  19. b=[50000;2250;-500;-2000];
  20. lb=zeros(3,1);
  21. [x,fval]=intlinprog(c,intcon,A,b,[],[],lb)
  23. %% 例3
  24. c=[-3;-2;-1]; intcon=3; % x3限定为整数
  25. A=ones(1,3); b=7;
  26. Aeq=[4 2 1]; beq=12;
  27. lb=zeros(3,1); ub=[+inf;+inf;1]; %x(3)为0-1变量
  28. [x,fval]=intlinprog(c,intcon,A,b,Aeq,beq,lb,ub)

（2）Matlab中线性 0 - 1 规划求解

  • 0 -1 规划无非就是把上界定为1，下届定为0的整数规划 （限制 lb、un）

在这里插入图片描述

七、整数规划的典型例题

（1）例题一：背包问题

在这里插入图片描述

  1. %% 背包问题（货车运送货物的问题）
  2. c = -[540 200 180 350 60 150 280 450 320 120];  % 目标函数的系数矩阵(最大化问题记得加负号)
  3. intcon=[1:10];  % 整数变量的位置(一共10个决策变量，均为0-1整数变量)
  4. A = [6 3 4 5 1 2 3 5 4 2];  b = 30;   % 线性不等式约束的系数矩阵和常数项向量（物品的重量不能超过30）
  5. Aeq = []; beq =[];  % 不存在线性等式约束
  6. lb = zeros(10,1);  % 约束变量的范围下限
  7. ub = ones(10,1);  % 约束变量的范围上限
  8. %最后调用intlinprog()函数
  9. [x,fval]=intlinprog(c,intcon,A,b,Aeq,beq,lb,ub)
  10. fval = -fval

（2）例题二：指派问题

在这里插入图片描述

  1. %% 指派问题（选择队员去进行游泳接力比赛）
  2. clear;clc
  3. c = [66.8 75.6 87 58.6 57.2 66 66.4 53 78 67.8 84.6 59.4 70 74.2 69.6 57.2 67.4 71 83.8 62.4]';  % 目标函数的系数矩阵（先列后行的写法）
  4. intcon = [1:20];  % 整数变量的位置(一共20个决策变量，均为0-1整数变量)
  5. % 线性不等式约束的系数矩阵和常数项向量（每个人只能入选四种泳姿之一，一共五个约束）
  6. A = [1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
  7.        0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0;
  8.        0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0;
  9.        0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0;
  10.        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1];
  11. % A = zeros(5,20);
  12. % for i = 1:5
  13. %     A(i, (4*i-3): 4*i) = 1;
  14. % end
  15. b = [1;1;1;1;1];
  16. % 线性等式约束的系数矩阵和常数项向量 （每种泳姿有且仅有一人参加，一共四个约束）
  17. Aeq = [1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0;
  18.           0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0;
  19.           0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0;
  20.           0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1];
  21. % Aeq = [eye(4),eye(4),eye(4),eye(4),eye(4)];  % 或者写成 repmat(eye(4),1,5)  
  22. beq = [1;1;1;1];
  23. lb = zeros(20,1);  % 约束变量的范围下限
  24. ub = ones(20,1);  % 约束变量的范围上限
  25. %最后调用intlinprog()函数
  26. [x,fval] = intlinprog(c,intcon,A,b,Aeq,beq,lb,ub)
  27. % reshape(x,4,5)' 矩阵变形
  28. %      0     0     0     1    甲自由泳
  29. %      1     0     0     0    乙蝶泳
  30. %      0     1     0     0    丙仰泳
  31. %      0     0     1     0    丁蛙泳
  32. %      0     0     0     0    戊不参加

（3）例题三：钢管切割问题

在这里插入图片描述

  1. %% 钢管切割问题
  2. %% (1)枚举法找出同一个原材料上所有的切割方法
  3. for i = 0: 2  % 2.9m长的圆钢的数量
  4.     for j = 0: 3  % 2.1m长的圆钢的数量
  5.         for k = 0:6   % 1m长的圆钢的数量
  6.             if 2.9*i+2.1*j+1*k >= 6 && 2.9*i+2.1*j+1*k <= 6.9
  7.                 disp([i, j, k])
  8.             end
  9.         end
  10.     end
  11. end
  12. % 老的MATLAB版本，会出现浮点数计算的误差
  13. % 只需要将上面的if这一行进行适当的放缩即可。
  14. % if 2.9*i+2.1*j+1*k >= 6-0.0000001 && 2.9*i+2.1*j+1*k <= 6.9+0.0000001
  16. %% (2) 线性整数规划问题的求解
  17. c = ones(7,1);  % 目标函数的系数矩阵
  18. intcon=[1:7];  %  整数变量的位置(一共7个决策变量，均为整数变量)
  19. A = -[1 2 0 0 0 0 1;  
  20.          0 0 3 2 1 0 1;
  21.          4 1 0 2 4 6 1];  % 线性不等式约束的系数矩阵
  22. b = -[100 100 100]'; %  线性不等式约束的常数项向量
  23. lb = zeros(7,1); % 约束变量的范围下限
  24. [x,fval]=intlinprog(c,intcon,A,b,[],[],lb)

八、最大化最小化模型

（1）模型的一般形式

在这里插入图片描述

  1. [x,feval] = fminimax(@Fun,x0,A,b,Aeq,beq,lb,ub，@nonlfun,option)

（2）典型例题

在这里插入图片描述

（3）模型的求解

  1. [x,feval] = fminimax(@Fun,x0,A,b,Aeq,beq,lb,ub，@nonlfun,option)

在这里插入图片描述

  • 主函数

    %% 最大最小化模型 : min{max[f1,f2,···,fm]}
    x0 = [6, 6]; % 给定初始值
    lb = [3, 4]; % 决策变量的下界
    ub = [8, 10]; % 决策变量的上界
    [x,feval] = fminimax(@Fun,x0,[],[],[],[],lb,ub)
    max(feval)
    % x =
    % 8.0000 8.5000
    % feval =
    % 13.5000 5.5000 5.5000 12.5000 8.5000 8.5000 5.5000 13.5000 9.5000 0.5000
    % 结论：
    % 在坐标为(8,8.5)处建立供应中心可以使该点到各需求点的最大距离最小，最小的最大距离为13.5单位。

  • 被调用的函数：Fun

    function f = Fun(x)

    1. a=[1 4 3 5 9 12 6 20 17 8];
    2. b=[2 10 8 18 1 4 5 10 8 9];
    3. %  函数向量
    4. f=zeros(10,1);
    5. for i = 1:10
    6.     f(i) = abs(x(1)-a(i))+abs(x(2)-b(i));  
    7. end

    % f(1) = abs(x(1)-a(1))+abs(x(2)-b(1));
    % f(2) = abs(x(1)-a(2))+abs(x(2)-b(2));
    % f(3) = abs(x(1)-a(3))+abs(x(2)-b(3));
    % f(4) = abs(x(1)-a(4))+abs(x(2)-b(4));
    % f(5) = abs(x(1)-a(5))+abs(x(2)-b(5));
    % f(6) = abs(x(1)-a(6))+abs(x(2)-b(6));
    % f(7) = abs(x(1)-a(7))+abs(x(2)-b(7));
    % f(8) = abs(x(1)-a(8))+abs(x(2)-b(8));
    % f(9) = abs(x(1)-a(9))+abs(x(2)-b(9));
    % f(10) = abs(x(1)-a(10))+abs(x(2)-b(10));
    end

九、多目标规划模型

（1）求解思路

  1. 先统一量纲和讲目标函数统一为最大化和最小化；
  2. 通过加权求和把多目标规划转换为单目标规划
  3. 求解得到符号约束条件的最优解，可以再通过最优解还原多目标函数
    在这里插入图片描述

（2）典型例题

在这里插入图片描述
在这里插入图片描述

  1. %%  多目标规划问题
  2. w1 = 0.4;  w2 = 0.6;  % 两个目标函数的权重  x1 = 5  x2 = 2
  3. w1 = 0.5;  w2 = 0.5;  % 两个目标函数的权重  x1 = 5  x2 = 2
  4. w1 = 0.3;  w2 = 0.7;  % 两个目标函数的权重  x1 = 1  x2 = 6
  5. c = [w1/30*2+w2/2*0.4 ;w1/30*5+w2/2*0.3];  % 线性规划目标函数的系数
  6. A = [-1 -1];  b = -7; % 不等式约束
  7. lb = [0 0]'; ub = [5 6]'; % 上下界
  8. [x,fval] = linprog(c,A,b,[],[],lb,ub)
  9. f1 = 2*x(1)+5*x(2)
  10. f2 = 0.4*x(1) + 0.3*x(2)
  13. %% 敏感性分析
  14. clear;clc
  15. W1 = 0.1:0.001:0.5;  W2 = 1- W1;  
  16. n =length(W1);
  17. F1 = zeros(n,1);  F2 = zeros(n,1);   X1 = zeros(n,1);  X2 = zeros(n,1);   FVAL = zeros(n,1);
  18. A = [-1 -1];  b = -7; % 不等式约束
  19. lb = [0 0]; ub = [5 6]; % 上下界
  20. for i = 1:n
  21.     w1 = W1(i);  w2 = W2(i);
  22.     c = [w1/30*2+w2/2*0.4 ;w1/30*5+w2/2*0.3];  % 线性规划目标函数的系数
  23.     [x,fval] = linprog(c,A,b,[],[],lb,ub);
  24.     F1(i) = 2*x(1)+5*x(2);
  25.     F2(i) = 0.4*x(1) + 0.3*x(2);
  26.     X1(i) = x(1);
  27.     X2(i) = x(2);
  28.     FVAL(i) = fval;
  29. end
  31. % 「Matlab」“LaTex字符汇总”讲解：https://blog.csdn.net/Robot_Starscream/article/details/89386748
  32. % 在图上可以加上数据游标，按住Alt加鼠标左键可以设置多个数据游标出来。
  33. figure(1) 
  34. plot(W1,F1,W1,F2)
  35. xlabel('f_{1}的权重') 
  36. ylabel('f_{1}和f_{2}的取值')
  37. legend('f_{1}','f_{2}')
  39. figure(2)
  40. plot(W1,X1,W1,X2)
  41. xlabel('f_{1}的权重') 
  42. ylabel('x_{1}和x_{2}的取值')
  43. legend('x_{1}','x_{2}')
  45. figure(3)
  46. plot(W1,FVAL)  % 看起来是两个直线组合起来的下半部分
  47. xlabel('f_{1}的权重') 
  48. ylabel('综合指标的值')
  52. % % 注意：代码文件仅供参考，一定不要直接用于自己的数模论文中
  53. % % 国赛对于论文的查重要求非常严格，代码雷同也算作抄袭

十、拓展例题

（1）护士值班问题

在这里插入图片描述
在这里插入图片描述

  1. %% 护士排班问题
  2. clear 
  3. % (1) 系数向量
  4. c = ones(6,1); 
  5. % (2) 不等式约束
  6. A =zeros(6,6);
  7. A(1,1) = -1; A(1,6) = -1;
  8. for i = 1:5
  9.     A(i+1, i) = -1;  A(i+1,i+1) = -1;
  10. end
  11. b = -[60 70 60 50 20 30]';
  12. %（3）上下界
  13. lb = zeros(6,1);
  15. % 注意，这题应该是一个整数规划问题哦
  16. intcon = [1:6];
  17. [x fval] = intlinprog(c,intcon, A, b, [], [], lb)  
  18. % fval =
  19. %    150
  20. % 注：本题的x可能有多个解！

（2）货舱装运问题

在这里插入图片描述

在这里插入图片描述

  1. %% 货机装货问题
  2. clear
  3. format long g   %可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）
  4. % 输入已知数据
  5. w = [10 16 8]'; % 前、中、后仓的重量限制
  6. v = [6800 8700 5300]'; % 前、中、后仓的体积限制
  7. a = [18 15 23 12]';  % 货物的重量
  8. b = [480 650 580 390]';  % 货物单位重量占用的空间
  9. c = [3100 3800 3500 2850]'; %  货物单位重量的利润
  10. % (1) 系数向量
  11. cc = [];  % 因为c这个符号被占用了，因此我们用cc表示系数向量
  12. for i = 1:4
  13.     for j = 1:3
  14.         cc = [cc;c(i)];  % 用循环的方法生成cc
  15.     end
  16. end
  17. cc = -cc;  % 因为要求最大值，所以这里对系数向量变号
  18. % (2) 不等式约束
  19. A = zeros(10,12); % 10个不等式约束
  20. for i = 1:4
  21.     A(i,3*i-2) = 1; A(i,3*i-1) = 1;  A(i,3*i) = 1;
  22. end
  23. for i = 5:7
  24.      j = i-4;
  25.     A(i,j) = 1; A(i,j+3) = 1;  A(i,j+6) = 1; A(i,j+9) = 1;
  26. end
  27. for i = 8:10
  28.      j = i-7;
  29.     A(i,j) = b(1); A(i,j+3) = b(2);  A(i,j+6) = b(3); A(i,j+9) = b(4);
  30. end
  31. b = [a;w;v];
  32. % (3) 等式约束
  33. Aeq = zeros(2,12);
  34. for i = 1:4
  35.     Aeq(1,3*i-2) = 1/10;
  36.     Aeq(1,3*i-1) = -1/16;
  37. end
  38. for i = 1:4
  39.     Aeq(2,3*i-2) = 1/10;
  40.     Aeq(2,3*i) = -1/8;
  41. end
  42. beq = [0,0]';
  43. %（4）上下界
  44. lb = zeros(12,1);
  45. % 进行求解
  46. [x fval] = linprog(cc, A, b, Aeq, beq, lb)
  47. fval = -fval
  48. % fval =
  49. %           121515.789473684
  50. % 注：本题的x可能有多个解！
  52. % 下面将x重新排列成4*3的矩阵
  53. % 注意这里不能使用reshape函数，reshape函数是按照列的顺序进行转换的，也就是第一列读完，读第二列，即x1对应x_1,1，x2对应x_2,1
  54. % 我们这题之前在定义x1-x12时，x1对应的是x_1,1, x2对应的是x_1,2 
  55. % 我们可以这样操作：先变形后转置
  56. x = reshape(x,3,4)'
  58. % 如果不能理解上面的操作，可以看下面这个简单的例题
  59. % y = [1,2,3,4,5,6];  % 想要把y里面的元素按照行的顺序来转换成大小为3*2的矩阵
  60. % yy = reshape(y,2,3)  % 先按照列的顺序转换为2*3的矩阵
  61. % yyy = yy'  % 对上一步操作得到的结果来进行转置

（3）非线性规划的求解

在这里插入图片描述

  • 主函数

    %% 拓展三：非线性规划的求解
    clear;clc
    A=[-1 -2 0]; b=-1; % 线性不等式约束
    x0 = [1 0 1]; % 初始值
    lb = [0 -inf -inf]; % 下界
    [x,fval] = fmincon(@fun,x0,A,b,[],[],lb,[],@nonlfun)
    fval = -fval

    % 先使用蒙特卡罗算法确定初始值，然后再进行计算
    N = 1000000;
    x1 = unifrnd(0,3,N,1); % x1在0~3之间均匀分布
    x2 = unifrnd(-8,7,N,1); % x2在-8~7之间均匀分布
    x3 = 2-x1.^2; % x1^2 + x3 = 2, 注意这里要用“.^”，表示每个元素分别平方
    fmax=-inf; % 初始化函数f的最大值为负无穷（后续只要找到一个比它大的我们就对其更新）
    for i=1:N

    1. x = [x1(i), x2(i), x3(i)];  %构造x向量, 这里千万别写成了：x =[x1, x2, x3]
    2. if (x(1)+2*x(1)^2+x(2)+2*x(2)^2+x(3)-10<=0) & (x(1)+x(1)^2+x(2)+x(2)^2-x(3)-50<=0) & (2*x(1)+x(1)^2+2*x(2)+x(3)-40<=0) & (x(1)+2*x(2)>=1)   % 判断是否满足条件
    3.     result =2*x(1)+3*x(1)^2+3*x(2)+x(2)^2+x(3) ;  % 如果满足条件就计算函数值
    4.     if  result  > fmax  % 如果这个函数值大于我们之前计算出来的最大值
    5.         fmax = result;  % 那么就更新这个函数值为新的最大值
    6.         x0 = x;  % 并且将此时的x1 x2 x3更新为初始值
    7.     end
    8. end

    end
    disp(‘蒙特卡罗算法选取的初始值为：’); disp(x0)
    [x,fval] = fmincon(@fun,x0,A,b,[],[],lb,[],@nonlfun)
    fval = -fval

  • 被调用函数：fun （目标函数）

    function f=fun(x)

    1. f=2*x(1)+3*x(1)^2+3*x(2)+x(2)^2+x(3);
    2. f=-f;  % 最大值加负号改为最小值

    end

  • 被调用函数：nonlfun （非线性约束条件）

    function [c,ceq]=nonlfun(x)

    1. % 非线性不等式约束    
    2. c=[x(1)+2*x(1)^2+x(2)+2*x(2)^2+x(3)-10
    3.     x(1)+x(1)^2+x(2)+x(2)^2-x(3)-50
    4.     2*x(1)+x(1)^2+2*x(2)+x(3)-40];
    5. % 非线性等式约束    
    6. ceq=x(1)^2+x(3)-2;

    end

（4）覆盖问题

在这里插入图片描述

  1. %%  第四题：覆盖问题
  2. c = ones(1,6);  % 目标函数的系数矩阵
  3. intcon=[1:6];  % 整数变量的位置(一共6个决策变量，均为0-1整数变量)
  4. % 线性不等式约束的系数矩阵和常数项向量
  5. A =- [1 1 1 0 0 0;
  6.          0 1 0 1 0 0;
  7.          0 0 1 0 1 0;
  8.          0 0 0 1 0 1;
  9.          0 0 0 0 1 1;
  10.          1 0 0 0 0 0;
  11.          0 1 0 1 0 1];
  12. b = -ones(7,1);
  13. Aeq = []; beq = [];  % 不存在线性等式约束
  14. lb = zeros(1,6);  % 约束变量的范围下限
  15. ub = ones(1,6);  % 约束变量的范围上限
  16. %最后调用intlinprog()函数
  17. [x,fval]=intlinprog(c,intcon,A,b,Aeq,beq,lb,ub)

发表评论

表情：
评论列表 （有 0 条评论，18人围观）

还没有评论，来说两句吧...

相关阅读