基本上这是一个开锁问题,锁的销有 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
——我不知道你的目标函数的细节,所以可能会有更多。例如,我目前使用 simpleround()
将输入转换X
为整数,但使用'PopulationType', 'custom'
(使用 custom'CreationFcn'
等'MutationFcn'
)可能会产生更好的结果。此外,'Vectorized'
可能会加快速度,但我不知道您的函数是否易于矢量化。
是的,我使用嵌套函数(我就是喜欢这些东西!);如果您使用子函数或独立函数,它会阻止这些巨大的、通常相同的输入参数列表,并且它们确实可以提高性能,因为几乎没有数据复制。但是,我意识到它们的作用域规则使它们有点类似于goto
构造,因此它们是-嗯-“不是每个人的一杯茶”……您可能希望将它们转换为子功能,以防止与您的冗长而无用的讨论合作伙伴 :)
无论如何,这应该是一个很好的起点。让我知道这是否有用。