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

基本上这是一个开锁问题，锁的销有 20 个不同的位置，有 12 个销。还有：

• 某些引脚的位置将被阻止，具体取决于所有其他引脚的位置。
• 根据锁具的具体情况，可能有多个适合的钥匙

...有趣！

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

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

可能的向量s（给予或接受一些，没有考虑过+1等）。对您相对简单的f(s) 进行约 7.4 亿次评估不应该超过 2 秒，并且在找到最大化f(s) 的所有s后，您将面临寻找线性组合的问题在您的向量集中，加起来就是其中一个解决方案s。

当然，这种组合的发现绝非易事，如果你正在处理，整个方法无论如何都会崩溃

int16:   ans = 2.311325368800510e+018
int32:   ans = 4.253529737045237e+037
int64:   ans = 1.447401115466452e+076

所以，我将在这里讨论一种更直接、更通用的方法。

由于我们讨论的是整数和相当大的搜索空间，我建议使用分支定界算法。但与bintprog 算法不同，您必须使用不同的分支策略，当然，这些策略应该基于非线性目标函数。

不幸的是，优化工具箱（或据我所知的文件交换）中没有这样的东西。 fmincon 是不行的，因为它使用梯度和 Hessian 信息（整数通常全为零），而 fminsearch 是不行的，因为你真的需要一个 em> 良好的初始估计，收敛速度（大致）为O(N)，这意味着对于这个 20 维问题，你必须等待很长时间才能收敛，没有保证找到了全局解决方案。

interval method 可能是可能的，但是，我个人对此没有什么经验。在 MATLAB 或其任何工具箱中没有与本机区间相关的东西，但有免费可用的 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)。 p>

特别注意，这目前（可能）是使用ga 的次优方式——我不知道你的目标函数的细节，所以可能会有更多。例如，我目前使用简单的round() 将输入X 转换为整数，但使用'PopulationType', 'custom'（带有自定义'CreationFcn'、'MutationFcn' 等）可能会产生更好的结果。此外，'Vectorized' 可能会加快速度，但我不知道您的函数是否易于矢量化。

是的，我使用嵌套函数（我就是喜欢这些东西！）；如果您使用子函数或独立函数，它会阻止这些巨大的、通常相同的输入参数列表，并且它们确实可以提高性能，因为几乎没有数据复制。但是，我意识到它们的作用域规则使它们有点类似于goto 构造，所以它们是-ahum-“不是每个人的一杯茶”......您可能希望将它们转换为子函数以防止冗长和无用与您的同事讨论:)

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

Answer 2

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

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

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

Answer 3

我知道，这个答案真的很适合您late。

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

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

这样，可以引导蛮力。