如何使用匿名功能优化MATLAB中的约束积分表达式？

Question

我有一个综合错误表达式 E = int[ abs(x-p)^2 ]dx 限制 x|0 到 x| L。变量p 是一个多项式，形式为 2*(a*sin(x)+b(a)*sin(2*x)+c(a)*sin(3*x))时间>。换句话说，系数b和c都是a的已知表达式。一个额外的方程被给出为 dE/da = 0。如果定义了上限L，则方程组是封闭的，我可以求解a，给出三个系数。

我设法获得了一个优化例程来解决 a 纯粹基于最大化 L。这可以通过在下面的代码中设置optimize=0 来确认。它提供的解决方案与我分析解决问题的解决方案相同。因此，我知道求解系数a 的方程是正确的。

我知道我提出的示例可以用铅笔和纸解决，但我正在尝试构建一个针对此类问题通用的优化函数（我有很多要评估的内容）。理想情况下，polynomial 作为函数的输入参数给出，然后输出xsol。显然，在我担心泛化之前，我需要对我在此处介绍的polynomial 进行优化。

无论如何，我现在需要通过一些约束进一步优化问题。首先，选择L。这让我可以计算a。一旦知道a，polynomial 就只是x 的已知函数，即p(x)。然后我需要从0->x 中确定满足以下约束的最大INTERVAL：|dp(x)/dx - 1| < tol。这让我可以用系数a 衡量polynomial 的性能。间隔就是我所说的“带宽”。我想强调两点：1）“带宽”与L 不同。 2）“带宽”内x的所有值必须满足约束。函数dp(x)/dx 确实会在公差标准内外波动，因此测试x 单个值的标准不起作用。它必须在一段时间内进行测试。第一个违规实例定义了带宽。我需要最大化这个“带宽”/间隔。对于输出，我还需要知道哪个L 导致了这种优化，因此我知道为给定约束选择正确的a。那是正式的问题陈述。 （希望这次我做对了）

现在我的问题是使用 MATLAB 的优化工具来设置这一切。我尝试遵循以下文章中的想法：

• Tutorial for the Optimization Toolbox™

为if statement 设置optimize=1 将与约束优化一起使用。我认为一些嵌套优化是如何涉及的，但我无法得到任何工作。我从 IMSL 优化库中提供了该问题的已知解决方案以进行比较/检查。它们写在优化例程下面。无论如何，这是我迄今为止整理的代码：

function [history] = testing()

% History
history.fval = [];
history.x = [];
history.a = []; 

%----------------
% Equations
polynomial = @(x,a) 2*sin(x)*a + 2*sin(2*x)*(9/20 -(4*a)/5) + 2*sin(3*x)*(a/5 - 2/15);
dpdx = @(x,a) 2*cos(x)*a + 4*cos(2*x)*(9/20 -(4*a)/5) + 6*cos(3*x)*(a/5 - 2/15);

% Upper limit of integration
IC = 0.8;       % initial
LB = 0;         % lower
UB = pi/2;      % upper

% Optimization
tol = 0.003;


% Coefficient
% --------------------------------------------------------------------------------------------
dpda = @(x,a) 2*sin(x) + 2*sin(2*x)*(-4/5) + 2*sin(3*x)*1/5;
dEda = @(L,a) -2*integral(@(x) (x-polynomial(x,a)).*dpda(x,a),0,L);
a_of_L = @(L) fzero(@(a)dEda(L,a),0);                                   % Calculate the value of "a" for a given "L"
EXITFLAG = @(L) get_outputs(@()a_of_L(L),3);                            % Be sure a zero is actually calculated


% NL Constraints
% --------------------------------------------------------------------------------------------
% Equality constraint (No inequality constraints for parent optimization)
ceq = @(L) EXITFLAG(L) - 1; % Just make sure fzero finds unique solution 
confun = @(L) deal([],ceq(L));

% Objective function
% --------------------------------------------------------------------------------------------
% (Set optimize=0 to test coefficent equations and proper maximization of L )
optimize = 1;

if optimize

%%%%  Plug in solution below

else 
    % Optimization options
    options = optimset('Algorithm','interior-point','Display','iter','MaxIter',500,'OutputFcn',@outfun);

    % Optimize objective
    objective = @(L) -L;
    xsol = fmincon(objective,IC,[],[],[],[],LB,UB,confun,options);

    % Known optimized solution from IMSL library
    % a = 0.799266;
    % lim = pi/2;
    disp(['IMSL coeff (a): 0.799266     Upper bound (L): ',num2str(pi/2)])
    disp(['code coeff (a): ',num2str(history.a(end)),'   Upper bound: ',num2str(xsol)])

end




    % http://stackoverflow.com/questions/7921133/anonymous-functions-calling-functions-with-multiple-output-forms
    function varargout = get_outputs(fn, ixsOutputs)
        output_cell = cell(1,max(ixsOutputs));
        [output_cell{:}] = (fn());
        varargout = output_cell(ixsOutputs);
    end

    function stop = outfun(x,optimValues,state)
        stop = false;

        switch state
            case 'init'
            case 'iter'
                % Concatenate current point and objective function
                % value with history. x must be a row vector.
                history.fval = [history.fval; optimValues.fval];
                history.x = [history.x; x(1)];
                history.a = [history.a; a_of_L(x(1))];

            case 'done'
            otherwise
        end
    end
end

我真的可以使用一些帮助来设置约束优化。我不仅对优化不熟悉，而且从未使用过 MATLAB。我还应该注意，我上面的内容不起作用，并且对于约束优化是不正确的。

更新：我在 if optimize 部分添加了一个 for 循环，以显示我试图通过优化实现的目标。显然，我可以只使用它，但它似乎非常低效，特别是如果我增加range 的分辨率并且必须多次运行此优化。如果您取消注释这些图，它将显示bandwidth 的行为方式。通过在整个范围内循环，我基本上测试了每个L，但肯定有更有效的方法来做到这一点？？

更新：已解决

Answer 1

所以看来fmincon 并不是完成这项工作的唯一工具。事实上，我什至无法让它工作。下面，fmincon 被“卡”在 IC 上，拒绝做任何事情……为什么……那是另一个帖子！使用相同的布局和公式，fminbnd 找到了正确的解决方案。据我所知，唯一的区别是前者使用了条件。但我的条件并不花哨，而且真的不需要。所以它必须与算法有关。我想这就是使用“黑匣子”时得到的结果。无论如何，经过漫长而痛苦的学习经历，这里有一个解决方案：

options = optimset('Display','iter','MaxIter',500,'OutputFcn',@outfun);

% Conditional
index = @(L) min(find(abs([dpdx(range(range<=L),a_of_L(L)),inf] - 1) - tol > 0,1,'first'),length(range));

% Optimize
%xsol = fmincon(@(L) -range(index(L)),IC,[],[],[],[],LB,UB,confun,options);
xsol = fminbnd(@(L) -range(index(L)),LB,UB,options);

我要特别感谢@AndrasDeak 的所有支持。如果没有帮助，我不会弄明白的！