Joh*_*nes 9 matlab mathematical-optimization
我有12组向量(每组约10-20个向量),我想选择每组的一个向量,以便将这些向量的总和作为参数的函数f最大化.另外,我对该总和的某些组件有约束.
例:
a_1 = [3 2 0 5], a_2 = [3 0 0 2], a_3 = [6 0 1 1], ... , a_20 = [2 12 4 3]
b_1 = [4 0 4 -2], b_2 = [0 0 1 0], b_3 = [2 0 0 4], ... , b_16 = [0 9 2 3]
...
l_1 = [4 0 2 0], l_2 = [0 1 -2 0], l_3 = [4 4 0 1], ... , l_19 = [3 0 9 0]
s = [s_1 s_2 s_3 s_4] = a_x + b_y + ... + l_z
约束:
s_1 > 40
s_2 < 100
s_4 > -20
目标:选择x,y,...,z以最大化f(s):
f(s) -> max
其中f是一个非线性函数,它取矢量s并返回一个标量.
强制执行需要太长时间,因为有大约5.9万亿个组合,并且由于我需要最大(甚至更好的前10个组合),我不能使用任何我认为的贪婪算法.
向量非常稀疏,大约70-90%是零.如果这有助于某种方式...？
Matlab优化工具箱没有任何帮助,因为它不支持离散优化.
基本上这是一个锁定拾取问题,锁的引脚有20个不同的位置,有12个引脚.也:
...有趣!
基于Rasman的方法和Phpdna的评论,以及假设您使用的int8数据类型,在给定的约束条件下有
>> d = double(intmax('int8'));
>> (d-40) * (d+100) * (d+20) * 2*d
ans =
737388162
可能的向量s(给予或采取一些,没有考虑+ 1等).相对简单的约7.4亿次评估f(s)不应超过2秒,并且找到所有 s最大化后f(s),您将面临在向量集中找到线性组合的问题,这些组合可以添加到其中一个解决方案中s.
当然,这种组合的发现并非易事,如果你正在处理,整个方法无论如何都会崩溃
int16: ans = 2.311325368800510e+018
int32: ans = 4.253529737045237e+037
int64: ans = 1.447401115466452e+076
所以,我将在这里讨论一种更直接,更通用的方法.
由于我们在谈论整数和相当大的搜索空间,我建议使用分支定界算法.但与bintprog算法不同,您必须使用不同的分支策略,当然,这些应该基于非线性目标函数.
不幸的是,在优化工具箱(或我可以找到的文件交换)中没有这样的东西.fmincon是不行的,因为它使用梯度和Hessian信息(对于整数通常都是零),并且fminsearch是不可行的,因为你需要一个非常好的初始估计,并且收敛速度是(大概)O(N),意思是,对于这个20维问题,你必须在收敛之前等待很长时间,而不保证找到全局解决方案.
一个间隔方法可能是一种可能性,但是,我个人有这样的经验非常少.在MATLAB或其任何工具箱中没有与本地区间相关的东西,但是有免费的INTLAB.
因此,如果您不想实现自己的非线性二进制整数编程算法,或者不想用INTLAB进行冒险,那么实际上只剩下一件事:启发式方法.在此链接中存在类似的情况,并概述了解决方案:使用ga全局优化工具箱中的遗传算法().
我会像这样大致实现这个问题:
function [sol, fval, exitflag] = bintprog_nonlinear()
%// insert your data here
%// Any sparsity you may have here will only make this more
%// *memory* efficient, not *computationally*
data = [...
... %// this will be an array with size 4-by-20-by-12
... %// (or some permutation of that you find more intuitive)
];
%// offsets into the 3D array to facilitate indexing a bit
offsets = bsxfun(@plus, ...
repmat(1:size(data,1), size(data,3),1), ...
(0:size(data,3)-1)' * size(data,1)*size(data,2)); %//'
%// your objective function
function val = obj(X)
%// limit "X" to integers in [1 20]
X = min(max(round(X),1),size(data,3));
%// "X" will be a collection of 12 integers between 0 and 20, which are
%// indices into the data matrix
%// form "s" from "X"
s = sum(bsxfun(@plus, offsets, X*size(data,1) - size(data,1)));
%// XxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxX
%// Compute the NEGATIVE VALUE of your function here
%// XxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxXxX
end
%// your "non-linear" constraint function
function [C, Ceq] = nonlcon(X)
%// limit "X" to integers in [1 20]
X = min(max(round(X),1),size(data,3));
%// form "s" from "X"
s = sum(bsxfun(@plus, offsets, X(:)*size(data,1) - size(data,1)));
%// we have no equality constraints
Ceq = [];
%// Compute inequality constraints
%// NOTE: solver is trying to solve C <= 0, so:
C = [...
40 - s(1)
s(2) - 100
-20 - s(4)
];
end
%// useful GA options
options = gaoptimset(...
'UseParallel', 'always'...
...
);
%// The rest really depends on the specifics of the problem.
%// Useful to look at will be at least 'TolCon', 'Vectorized', and of course,
%// 'PopulationType', 'Generations', etc.
%// THE OPTIMZIATION
[sol, fval, exitflag] = ga(...
@obj, size(data,3), ... %// objective function, taking a vector of 20 values
[],[], [],[], ... %// no linear (in)equality constraints
1,size(data,2), ... %// lower and upper limits
@nonlcon, options); %// your "nonlinear" constraints
end
请注意,即使您的约束基本上是线性的,您必须计算值的方式也s需要使用自定义约束函数(nonlcon).
特别要注意的是,目前(可能)这是一种次优的使用方式ga- 我不知道你的目标函数的细节,所以可能还有更多.举例来说,我目前使用一个简单round()的输入转换X为整数,但使用'PopulationType', 'custom'(使用自定义'CreationFcn','MutationFcn'等)可能会产生更好的效果.此外,'Vectorized'可能会加快速度,但我不知道你的功能是否很容易被矢量化.
是的,我使用嵌套函数(我只喜欢那些东西!); 如果使用子函数或独立函数,它会阻止这些巨大的,通常相同的输入参数列表,并且它们实际上可以提高性能,因为几乎没有数据复制.但是,我意识到他们的范围规则使他们有点类似于goto结构,所以他们是-ahum-"不是每个人的一杯茶"...你可能想要将它们转换为子功能,以防止与你的长期和无用的讨论合作伙伴 :)
无论如何,这应该是一个很好的起点.让我知道这是否有用.