创建有效的优化问题

当线性问题具有整数约束时，solve 调用 intlinprog 来获得解。有关获得更快解或更多整数可行点的建议，请参阅调整整数线性规划。

在开始求解问题之前，有时您可以改进问题约束或目标的表示。通常，软件可以以向量化方式而不是循环方式更快地创建目标函数或约束的表达式。当优化表达式需要自动微分时，这种速度差异尤其大；请参阅 Optimization Toolbox 中的自动微分。

通常，将其写为单独的函数更容易创建高效循环，如创建 for 循环进行静态分析中所述。在这种情况下，使用 fcn2optimexpr 创建包含函数的优化表达式，如优化表达式的静态分析中所述。

假设您的目标函数是

$\sum_{i = 1}^{30} \sum_{j = 1}^{30} \sum_{k = 1}^{10} x_{i, j, k} b_{k} c_{i, j},$

其中 x、b 和 c 是优化变量。构造该目标函数的两种常用方法如下：

使用 for 循环。

x = optimvar('x',30,30,10);
b = optimvar('b',10);
c = optimvar('c',30,30);
tic
expr = optimexpr;
for i = 1:30
    for j = 1:30
        for k = 1:10
            expr = expr + x(i,j,k)*b(k)*c(i,j);
        end
    end
end
toc

Elapsed time is 307.459465 seconds.

这里，expr 包含目标函数表达式。为了让软件有效地评估这个表达式，将其写为一个单独的函数：

function expr = loopme(x,b,c)
expr = 0;
ni = size(x,1);
nj = size(x,2);
nk = size(x,3);
for i = 1:ni
    for j = 1:nj
        for k = 1:nk
            expr = expr + x(i,j,k)*b(k)*c(i,j);
        end
    end
end
end

包含使用 fcn2optimexpr 的表达式。

tic
expr = fcn2optimexpr(@loopme,x,b,c,OutputSize=[1,1]);
toc

Elapsed time is 23.888763 seconds.

在问题中包含目标函数。

problem = optimproblem("Objective",expr);

使用向量化语句。矢量化语句通常比未经静态分析修改的 for 循环运行得更快。您可以通过多种方式创建向量化语句。

展开 b 和 c。为了实现逐项乘法，请创建与 x 大小相同的常量。

tic
bigb = reshape(b,1,1,10);
bigb = repmat(bigb,30,30,1);
bigc = repmat(c,1,1,10);
expr = sum(sum(sum(x.*bigb.*bigc)));
toc

Elapsed time is 0.013631 seconds.

对 b 循环一次。

tic
expr = optimexpr;
for k = 1:10
    expr = expr + sum(sum(x(:,:,k).*c))*b(k);
end
toc

Elapsed time is 0.044985 seconds.

通过循环 b 然后在循环后对项求和来创建一个表达式。

tic
expr = optimexpr(30,30,10);
for k = 1:10
    expr(:,:,k) = x(:,:,k).*c*b(k);
end
expr = sum(expr(:));
toc

Elapsed time is 0.039518 seconds.

观察示例使用优化变量的约束静电非线性优化的向量化和非矢量化实现之间的速度差异。此示例使用 R2020b 中的自动微分进行计时。

N = 30;
x = optimvar('x',N,'LowerBound',-1,'UpperBound',1);
y = optimvar('y',N,'LowerBound',-1,'UpperBound',1);
z = optimvar('z',N,'LowerBound',-2,'UpperBound',0);
elecprob = optimproblem;
elecprob.Constraints.spherec = (x.^2 + y.^2 + (z+1).^2) <= 1;
elecprob.Constraints.plane1 = z <= -x-y;
elecprob.Constraints.plane2 = z <= -x+y;
elecprob.Constraints.plane3 = z <= x-y;
elecprob.Constraints.plane4 = z <= x+y;

rng default % For reproducibility
x0 = randn(N,3);
for ii=1:N
    x0(ii,:) = x0(ii,:)/norm(x0(ii,:))/2;
    x0(ii,3) = x0(ii,3) - 1;
end
init.x = x0(:,1);
init.y = x0(:,2);
init.z = x0(:,3);
opts = optimoptions('fmincon','Display','off');

tic
energy = optimexpr(1);
for ii = 1:(N-1)
    jj = (ii+1):N; % Vectorized
    tempe = (x(ii) - x(jj)).^2 + (y(ii) - y(jj)).^2 + (z(ii) - z(jj)).^2;
    energy = energy + sum(tempe.^(-1/2));
end
elecprob.Objective = energy;
disp('Vectorized computation time:')
[sol,fval,exitflag,output] = solve(elecprob,init,'Options',opts);
toc

Vectorized computation time:
Elapsed time is 1.838136 seconds.

tic
energy2 = optimexpr(1); % For nonvectorized comparison
for ii = 1:(N-1)
    for jjj = (ii+1):N; % Not vectorized
        energy2 = energy2 + ((x(ii) - x(jjj))^2 + (y(ii) - y(jjj))^2 + (z(ii) - z(jjj))^2)^(-1/2);
    end
end
elecprob.Objective = energy2;
disp('Non-vectorized computation time:')
[sol,fval,exitflag,output] = solve(elecprob,init,'Options',opts);
toc

Non-vectorized computation time:
Elapsed time is 204.615210 seconds.

向量化版本比非矢量化版本快约 100 倍。

与静态分析版本进行比较。myenergy 函数的代码在此示例末尾。

elecprob.Objective = fcn2optimexpr(@myenergy,N,x,y,z);

通过调用 solve 求解问题。对解进行计时。

disp('Static analysis computation time:')
tic
[sol,fval,exitflag,output] = solve(elecprob,init,'Options',opts);
toc

Static analysis computation time:
Elapsed time is 1.293513 seconds.

静态分析提供最少的计算时间。

使用以下代码创建 myenergy 函数。

function energy = myenergy(N, x, y, z)
energy = 0;
for ii = 1:(N-1)
    for jj = (ii+1):N
        energy = energy + ((x(ii) - x(jj))^2 + (y(ii) - y(jj))^2 + (z(ii) - z(jj))^2)^(-1/2);
    end
end
end

创建有效的优化问题

相关主题

WeChat