matlab - 在Matlab中使用离散参数进行优化

Question

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

Answer 1

基本上，这是一个 pry 锁问题，其中锁的销钉有20个不同的位置，并且有12个销钉。还:

该引脚的某些位置将被阻止，具体取决于所有其他引脚的位置。

根据锁的具体情况，可能有多个适合

的 key

...有趣的!

根据Rasman的方法和Phpdna的评论，并假设您使用int8作为数据类型，在给定的约束下，存在

>> d = double(intmax('int8'));
>> (d-40) * (d+100) * (d+20) * 2*d
ans =
    737388162

可能的向量s(给出或取几个，没有考虑+1等)。对您相对简单的f(s)进行的约7.4亿次评估不应花费超过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维问题，您必须等待很长时间才能收敛，而不能保证找到全局解决方案。

interval method是可能的，但是，我个人对此几乎没有经验。 MATLAB或其任何工具箱中都没有与 native 间隔相关的东西，但是有免费提供的INTLAB。

因此，如果您不想实现自己的非线性二进制整数编程算法，或者不想在INTLAB上冒险，那么实际上只剩下一件事了:启发式方法。在this link中，也有类似的情况，并给出了解决方案的概述:使用“全局优化”工具箱中的遗传算法(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-“不是每个人都可以喝到茶” ...您可能希望将它们转换为子功能，以防止与他们进行冗长而无用的讨论您的同事:)

无论如何，这应该是一个不错的起点。让我知道这是否有用。