[我分享] 建模训练案例-8:R2021b 优化工具箱更新简介

[复制链接]
qibbxxt 发表于 2021-9-29 14:30:59
MATLAB2021b优化/全局优化功能更新介绍

MATLAB2021b最近已经发布,Release Note中的更新涵盖几十个工具箱在不同应用场景下或大或小的改动,本篇文章重点介绍在优化/全局优化工具箱函数中,几个颇为令人意动,也是令我们期待已久的更新。

我们专栏文章所有优化问题的MATLAB求解代码,基本都是围绕全新问题式建模编写框架构思的,从R2017b到目前最新发布的R2021b,其中技巧用法细节在不同版本的不同处理方式,体现在专栏每一个优化模型的多解思路中,这一次也不例外,对之前版本操作方式的更新,多数仍然体现在这个框架的补充完善上,因此本篇文章重点将提到如下内容:

  • 全局优化工具箱对问题式建模流程的正式官方支持。优化工具箱函数在R2017b就已经支持问题式模型的构造、求解,但全局优化工具箱却仍然只支持求解器建模方式,使得用户如果想把问题式建模构造模型用ga、particleswarm、simulannealbnd等全局寻优函数求解,还要用prob2struct等函数做一次转换,新版本下,由于全局搜索函数调用方式支持问题式建模流程,这个困扰基本上彻底消除。本文首先将以ga为例,介绍利用问题式建模与全局优化函数之间结合的代码实现方式;
  • ga函数支持线性等式约束和整数变量的共存。之前版本如果问题出现等式约束就不能有整数变量,反之亦然,新版R2021b关于这一条件有了一定的放松:可以允许线性等式约束和整数变量的共存;
  • 问题式建模对非线性整数规划求解的支持。R2021b版本之前,问题式建模如果输入模型判定为非线性整数规划,则出现模型不可解的错误提示,现在则首先会通过ga尝试求解;
  • 求解函数solve对lsqlin/lsqnonlin类型问题的解析支持。R2021b版本之前,模型出现非线性项时,都默认选择fmincon求解,现在对于平方和形式的目标函数,只要满足对应约束类型的限定要求,可以自动解析为lsqlin或lsqnonlin进行求解。

全局优化工具箱问题式建模支持ga求解问题式建模模型-1(Till R2021a)

之前在专栏文章数学建模代码实战训练-4:料场寻址问题局部/全局寻优代码方案探析当中,探讨了利用全局优化工具箱命令如:ga、particleswarm等求解指定区域内的最佳送料场位置,撰文其时使用的是R2021a版本,文中提到全局优化工具箱不支持问题式建模。

为能让全局优化工具箱函数也能使用问题式建模构造出的模型,我们构思了借助prob2struct将问题式模型参数重新对应求解器模式模型struct域字段的方案,对模型预先实施转换,并在求解之后,用自编函数paser_result把求解优化结果重新转换为和optimvar维度一致的形式。

方便比较起见,以ga为例,重述其求解步骤:

  • 构造问题式建模模型,模型中整合了基础数据、设计变量、目标函数以及约束条件等,代码如下:

  1. function prob = ProbConstruct()
  2. LocWork = [1.25 8.75 .5 5.75 3 7.25;
  3.     1.25 .75 4.75 5 6.5 7.75];
  4. [S,L] = bounds(LocWork,2);
  5. ConsumeCem = [3 5 4 7 6 11]';
  6. Pos0 = [5 1;2 7];
  7. c = optimvar('c',6,2, 'LowerBound',0,'UpperBound',max(ConsumeCem)); % 构造设计变量
  8. x = optimvar('x',2,1, 'LowerBound',S(1),'UpperBound',L(1)); % 构造设计变量
  9. y = optimvar('y',2,1, 'LowerBound',S(2),'UpperBound',L(2)); % 构造设计变量
  10. f1=fcn2optimexpr(@(x,y)dist([x,y], LocWork)',x,y);
  11. prob = optimproblem('Objective',sum(sum(f1.*c))); % 构造决策目标
  12. prob.Constraints.con1 = sum(c,2)==ConsumeCem; % 构造约束条件
  13. prob.Constraints.con2 = sum(c)'<=[20;20]; % 构造约束条件
  14. end
复制代码
  • 调用ga求解。
  1. clear;clc;close all
  2. prob = ProbConstruct();
  3. Prostr = prob2struct(prob);
  4. opts = optimoptions('ga','PlotFcn',@gaplotbestf,...
  5.     'HybridFcn', @fmincon,...
  6.     'UseParallel', true,...
  7.     'PopulationSize', 100);
  8. [x,fval] = ga(Prostr.objective, numel(Prostr.lb),...
  9.     full(Prostr.Aineq), Prostr.bineq,...
  10.     full(Prostr.Aeq), Prostr.beq,...
  11.     Prostr.lb,Prostr.ub,[],Prostr.intcon, opts);
复制代码

在ga命令调用前,第3行已经用prob2struct把模型信息转换为结构体变量Prostr,因此第8~10行ga函数的主要输入参数均来自该变量不同的域字段。

  • 编写结果维度转换函数paser_result。由于求解结果中的变量(无论有多少个),在调用ga时被转换为一个单独的一维向量,这种形式的可读性是很低的,因此还需要把这些变量重新转换为原来问题式建模的形式,于是自编函数可以解决这一问题,代码如下:

  1. function sol = paser_result(prob, result, fvl)
  2. sol = struct;
  3. idx = varindex(prob);
  4. names = fieldnames(prob.Variables);
  5. for i = names'
  6.     iddx = getfield(idx, i{1}); %#ok<*GFLD>
  7.     v = reshape(result(iddx), size(getfield(prob.Variables, i{1})));
  8.     sol = setfield(sol, i{:}, v);
  9. end
  10. sol = setfield(sol, 'obj', fvl); %#ok<*SFLD>
  11. end
复制代码
  • 调用paser_result转换结果维度。
  1. sol = paser_result(prob, x, fval)
  2. sol =
  3.       c: [6×2 double]
  4.       x: [2×1 double]
  5.       y: [2×1 double]
  6.     obj: 93.7409
复制代码
从上述步骤看出,在R2021a或之前版本,想让全局优化工具箱命令`ga`求解以问题式建模构造的模型,需要先通过`prob2struct`转换为求解器建模的结构体形式,解析其中各个域的字段数据之后,才能用`ga`求解。

ga求解问题式建模模型-2(Since R2021b)

自R2021b版本,全局优化工具箱正式支持问题式建模的模型求解,这意味着以上代码中:问题式模型→求解器模型(prob2struct)→ga→结果维度转换(paser_result)的流程,可以被大幅简化,中间转换过程以及求解后处理部分的维度转换过程都可以省去,代码如下:

  1. clc;clear;close all;
  2. prob = ProbConstruct()
  3. opts = optimoptions('ga','PlotFcn',@gaplotbestf,...
  4.      'HybridFcn', @fmincon,...
  5.      'UseParallel', true,...
  6.      'PopulationSize', 100);
  7. [sol,fval] = solve(prob,"Solver","ga","Options",opts)
复制代码

运行结果:

  1. Starting parallel pool (parpool) using the 'local' profile ...
  2. Connected to the parallel pool (number of workers: 12).

  3. Solving problem using ga.
  4. Optimization terminated: average change in the fitness value less than options.FunctionTolerance.
  5. sol =
  6.     c: [6×2 double]
  7.     x: [2×1 double]
  8.     y: [2×1 double]
  9. fval = 93.7409
复制代码
ga等式约束的支持增强
截止R2021a版本,ga函数对于同时存在整数变量和等式约束的问题,使用存在一定的限制,即:等式约束和整数变量在调用ga求解时,只能出现其中之一,这在调用格式中也做了限定和扼要的说明。因此早期版本如果期望能解决同时具有这两个特征的问题,要把等式约束条件转换成两个等效的不等式约束,例如就需要转换为如下的不等式形式:





在R2021b版本中,一定程度上放宽了条件,允许整数变量和等式约束条件同时存在,因此这类问题的模型可以省去上述的变换操作,例如下列问题:


如不限定整数变量,则可以通过如下代码求解:

  1. x = optimvar('x',1,5, 'LowerBound',-5,'UpperBound',5); % 设计变量

  2. prob = optimproblem('Objective',(x(1)-1)^2+(x(1)-x(2))^2+(x(2)-x(3))^3+(x(3)-x(4))^4+(x(4)-x(5))^5); % 决策目标

  3. prob.Constraints.con1 = x(1)+x(2)+x(3)==3*sqrt(2)+2; % 约束条件1
  4. prob.Constraints.con2 = x(2)-x(3)+x(4)==2*sqrt(2)-2;   % 约束条件2
  5. opts = optimoptions('ga','PlotFcn',@gaplotbestf,...
  6.      'HybridFcn', @fmincon,...
  7.      'PopulationSize', 100);
  8. [sol1,fval1] = solve(prob,struct('x',randi(10,size(x))))  % 调用fmincon
  9. Solving problem using fmincon.

  10. Local minimum found that satisfies the constraints.

  11. Optimization completed because the objective function is non-decreasing in
  12. feasible directions, to within the value of the optimality tolerance,
  13. and constraints are satisfied to within the value of the constraint tolerance.

  14. <stopping criteria details>
  15. sol1 =     x: [5.0000 3.5355 -2.2929 -5.0000 5.0000]
  16. fval1 = -9.9730e+04

  17. [sol2,fval2] = solve(prob,"Solver","ga","Options",opts)  % 调用ga
  18. Solving problem using ga.
  19. Optimization terminated: maximum number of generations exceeded.
  20. sol2 =     x: [5.0000 3.5355 -2.2929 -5.0000 5.0000]
  21. fval2 = -9.9730e+04
复制代码
如果限定为整数变量,则fmincon无法求解,MATLAB官方工具箱只能通过ga求解此类非线性整数规划问题:

  1. clc;clear;close all;
  2. x = optimvar('x',1,4, 'LowerBound',-5,'UpperBound',5); % 设计变量
  3. xInt = optimvar('xInt',1,'type','integer', 'LowerBound',-5,'UpperBound',5); % 设计变量

  4. prob = optimproblem('Objective',(x(1)-1)^2+(x(1)-xInt)^2+(xInt-x(2))^3+(x(2)-x(3)^2)^4+(x(3)-x(4))^5); % 决策目标

  5. prob.Constraints.con1 = x(1)+xInt+x(2)==3*sqrt(2)+2; % 约束条件1
  6. prob.Constraints.con2 = xInt-x(2)+x(3)==2*sqrt(2)-2;   % 约束条件2
  7. opts = optimoptions('ga','PlotFcn',@gaplotbestf,...
  8.      'PopulationSize', 300);
  9. [sol3,fval3] = solve(prob,"Solver","ga","Options",opts)
  10. Solving problem using ga.
  11. Optimization terminated: average change in the penalty fitness value less than options.FunctionTolerance
  12. and constraint violation is less than options.ConstraintTolerance.
  13. sol3 =
  14.        x: [0.2377 1.0050 -3.1666 4.9126]
  15.     xInt: 5
  16. fval3 = -2.7709e+04
复制代码
注:目前整数变量只能和线性等式约束条件共存,如果是非线性等式约束则仍然需要采用之前的方式,转换为一组不等式约束才能使用ga求解。

由于问题比较简单,可以通过枚举方式,让变量xInt从​遍历,借助fmincon计算,查看结果的准确性:
  1. clc;clear;close all;
  2. x = optimvar('x',1,4, 'LowerBound',-5,'UpperBound',5); % 设计变量
  3. options = optimoptions('fmincon',...
  4.     'Algorithm','interior-point',...
  5.     'EnableFeasibilityMode',true,...
  6.     "Display","none",...
  7.     'SubproblemAlgorithm','cg');
  8. for xInt = -5:5
  9.     prob = optimproblem('Objective',(x(1)-1)^2+(x(1)-xInt)^2+(xInt-x(2))^3+(x(2)-x(3)^2)^4+(x(3)-x(4))^5); % 决策目标
  10.     prob.Constraints.con1 = x(1)+xInt+x(2)==3*sqrt(2)+2; % 约束条件1
  11.     prob.Constraints.con2 = xInt-x(2)+x(3)==2*sqrt(2)-2;   % 约束条件2
  12.     [solT,fvalT] = solve(prob,struct('x',randi(10,size(x))),'options',options)
  13. end
复制代码
注意:上述程序选择内点法算法,采用了`R2021b`一个新求解参数`'SubproblemAlgorithm'`,当该参数值为`'cg'`时,可以处理一部分用`fmincon`不容易找到可行解的情形,当然,这个问题在本例中并不存在,因此采用内点法或序列二次规划算法(参数`sqp`都可以很容易找到可行解)。之所以针对内点罚函数法,增加一个可行解求解增强选项,主要是算法层面的:内点罚函数法对于初值有严格要求,其必须落在可行区域内,在模型规模较大、变量较多时,用户很难找到这组初值。推测这是增加该选项的核心原因之一。

运行结果:
  1. solT =     x: [5.1995 2.6453 5.1952 4.2478]
  2. fvalT = 3.5092e+05
  3. solT =     x: [5.1996 2.2375 5.1959 3.4941]
  4. fvalT = 3.7569e+05
  5. solT =     x: [5.1997 2.6266 5.1961 4.0796]
  6. fvalT = 3.5279e+05
  7. solT =     x: [5.0400 2.6760 5.0398 4.3710]
  8. fvalT = 2.6661e+05
  9. solT =     x: [5.0000 2.2426 4.0711 4.9494]
  10. fvalT = 4.2197e+04
  11. solT =     x: [5.0000 1.2426 2.0711 4.9999]
  12. fvalT = -90.2773
  13. solT =     x: [5.0000 0.2427 0.0711 4.9998]
  14. fvalT = -2.8762e+03
  15. solT =     x: [5.0000 -0.7573 -1.9289 5.0000]
  16. fvalT = -1.5522e+04
  17. solT =     x: [3.9385 -0.6959 -2.8674 5.0000]
  18. fvalT = -2.3756e+04
  19. solT =     x: [2.0985 0.1441 -3.0275 5.0000]
  20. fvalT = -2.6648e+04
  21. solT =     x: [0.2513 0.9913 -3.1802 5.0000]
  22. fvalT = -2.9616e+04
复制代码

比较ga和fmincon的求解结果,可以看出遗传算法准确算出了整数变量的结果应为5,求解结果-27709略高于fmincon的最优解-29616。

优化工具箱

比较R2021a,优化工具箱的变化相对较少,除了之前提到的fmincon新参数'SubproblemAlgorithm',还有两个主要更新:

  • 当目标函数呈现:“”形式的时候,问题式建模流程会将问题类型解析为lsqlin(约束条件为线性),或lsqnonlin(变量有界);
  • 问题式建模求解时,对存在整数约束的非线性整数规划问题,以前版本会提示不可解,现在会自动尝试将其归类为遗传算法ga求解的问题。

更新1:lsqlin/lsqnonlin问题的解析支持

对于如下问题:

问题式建模对于此类目标函数平方和形式,且具有线性约束条件的模型求解,会自动归类为lsqlin问题:
  1. clc;clear;close all;
  2. x = optimvar('x',1,3, 'LowerBound',-5,'UpperBound',5); % 设计变量

  3. prob = optimproblem('Objective',(x(1)-1)^2+(x(2)+5)^2+(x(3)-2*x(1))^2); % 决策目标

  4. prob.Constraints.con1 = x(1)+x(2)<=3;           % 线性约束条件1
  5. prob.Constraints.con2 = x(1)-2*x(2)+x(3)<=12;   % 线性约束条件2
  6. [solL,fvalL] = solve(prob)
复制代码
求解结果信息:
  1. Solving problem using lsqlin.

  2. Minimum found that satisfies the constraints.

  3. Optimization completed because the objective function is non-decreasing in
  4. feasible directions, to within the value of the optimality tolerance,
  5. and constraints are satisfied to within the value of the constraint tolerance.

  6. <stopping criteria details>

  7. solL =  
  8.     x: [0.7857 -4.8571 1.5000]
  9. fvalL = 0.0714
复制代码
于目标非线性、变量有界(无其他约束)的模型,问题式建模求解归类为lsqnonlin问题:

  1. clc;clear;close all;
  2. x = optimvar('x',1,3, 'LowerBound',-5,'UpperBound',5); % 设计变量
  3. prob = optimproblem('Objective',(x(1)-1)^2+(x(2)+5)^2+(x(3)-2*sin(x(1)))^2); % 决策目标
  4. [solL,fvalL] = solve(prob,struct('x',randi(10,size(x))))
复制代码
注意目标函数整体是平方和形式,但变量带有正弦形式,其求解结果信息:
  1. Solving problem using lsqnonlin.

  2. Local minimum found.

  3. Optimization completed because the size of the gradient is less than
  4. the value of the optimality tolerance.

  5. <stopping criteria details>

  6. solL =
  7.     x: [1 -4.9994 1.6829]
  8. fvalL = 3.7253e-07
复制代码

更新2:非线性整数规划的支持增强

对于如下这种非线性整数规划问题,R2021a或以前的版本,问题式建模会提示错误,不能直接求解:


在新版本R2021b中,这种问题默认归类为使用ga求解:

  1. clc;clear;close all;
  2. x = optimvar('x',1,3, 'LowerBound',-5,'UpperBound',5); % 设计变量
  3. y = optimvar('y',1,'Type','integer');

  4. prob = optimproblem('Objective',(x(1)+12*y*x(3))^2+(x(2)+5*y)^2+(y*x(3)-2*x(1))^2);
  5. prob.Constraints.con1 = x(1)+x(2)+y*x(1)*x(3)<=6*y;
  6. prob.Constraints.con2 = x(1)-2*x(2)+y*x(3)>=5;

  7. [solG,fvalG] = solve(prob)
复制代码
求解结果:
  1. Solving problem using ga.
  2. Optimization terminated: average change in the penalty fitness value less than options.FunctionTolerance
  3. and constraint violation is less than options.ConstraintTolerance.

  4. solG =
  5.     x: [0.0461 -3.5730 -0.0045]
  6.     y: 1
  7. fvalG = 2.0456
复制代码

总结

从上述更新来看,此次R2021b版本在优化/全局优化工具箱的调整和更新,看似项目很少,但四两拨千斤,其真实力度并不小,首先,对问题式建模求解范围的扩大支持,实现了局部优化/全局优化工具箱求解方式的内部整合与逻辑自恰,不再是“局部优化一套打法、全局优化工具箱另一套打法”的尴尬情况;其次,扩展的类型支持也更符合优化理论问题的分类。







1 条回复


ilovematlab 发表于 2021-10-29 17:12:49
大家有兴趣也可以免费注册参加12月14日的《2021小迈步第七课 | MATLAB 轻松求解优化问题》,通过一些具体实例来学习如何调用 MATLAB 中丰富的优化算法来求解不同类型优化问题。
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表