9

我有 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 优化工具箱也没有帮助,因为它对离散优化的支持不多。

4

3 回答 3

6

基本上这是一个开锁问题,锁的销有 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构造,因此它们是-嗯-“不是每个人的一杯茶”……您可能希望将它们转换为子功能,以防止与您的冗长而无用的讨论合作伙伴 :)

无论如何,这应该是一个很好的起点。让我知道这是否有用。

于 2013-06-26T12:15:39.757 回答
0

除非你定义了一些关于向量集如何组织的智能,否则除了纯粹的蛮力之外,没有智能的方法可以解决你的问题。

假设您发现 s st f(s) 是给定 s 的最大约束,您仍然需要弄清楚如何使用 12 个 4 元素向量(如果曾经有一个超定系统)构建 s,其中每个向量有 20 个可能的值. 稀疏性可能会有所帮助,尽管我不确定如何使具有四个元素的向量70-90%为零,并且稀疏性只有在向量的组织方式方面有一些尚待描述的方法时才有用

所以我不是说你不能解决问题,我是说你需要重新考虑问题是如何设置的。

于 2013-06-26T14:12:13.693 回答
0

我知道,这个答案真的传到你了late

不幸的是,问题,就目前而言,没有多少模式可以利用,除了蛮力-Branch&Bound,Master&Slave等-尝试主从方法-即首先解决函数连续非线性问题作为主问题,然后解决离散的选择作为奴隶可能会有所帮助,但是有尽可能多的组合,并且没有关于向量的更多信息,没有太多的工作空间。

但是基于给定的几乎无处不在的连续函数,基于和和乘法运算符及其逆运算符的组合,稀疏性是这里要利用的一个明显点。如果 70-90% 的向量为零,则解空间的几乎很大一部分将接近于零,或接近infinite. 因此,一个 80-20 的伪解很容易丢弃“零”组合,而只使用“无限”组合。

这样,可以引导蛮力。

于 2015-04-12T05:49:23.200 回答