如何使用匿名函数优化 MATLAB 中的约束积分表达式? [英] How do I optimize constrained integral expressions in MATLAB using anonymous functions?

问题描述

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

I have an integrated error expression E = int[ abs(x-p)^2 ]dx with limits x|0 to x|L. The variable p is a polynomial of the form 2*(a*sin(x)+b(a)*sin(2*x)+c(a)*sin(3*x)). In other words, both coefficients b and c are known expressions of a. An additional equation is given as dE/da = 0. If the upper limit L is defined, the system of equations is closed and I can solve for a, giving the three coefficients.

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

I managed to get an optimization routine to solve for a purely based on maximizing L. This is confirmed by setting optimize=0 in the code below. It gives the same solution as if I solved the problem analytically. Therefore, I know the equations to solve for the coefficent a are correct.

我知道我提供的示例可以用铅笔和纸解决，但我正在尝试构建一个针对此类问题推广的优化函数(我有很多需要评估).理想情况下，polynomial 作为输入参数提供给一个函数，然后输出 xsol.显然，在我担心泛化之前，我需要让优化适用于我在这里介绍的 polynomial.

I know the example I presented can be solved with pencil and paper, but I'm trying to build an optimization function that is generalized for this type of problem (I have a lot to evaluate). Ideally, polynomial is given as an input argument to a function which then outputs xsol. Obviously, I need to get the optimization to work for the polynomial I presented here before I can worry about generalizations.

无论如何，我现在需要通过一些约束进一步优化问题.首先，选择 L.这允许我计算 a.一旦知道a，polynomial 只是x 的已知函数，即p(x).然后我需要从 0->x 确定满足以下约束的最大间隔:|dp(x)/dx - 1|<代码.这给了我一个带有系数 a 的 多项式 性能的度量.间隔就是我所说的带宽".我想强调两件事:1)带宽"与 L 不同.2)带宽"内的所有x值必须满足约束条件.函数 dp(x)/dx 确实会在容差标准内外振荡，因此测试 x 的单个值的标准不起作用.它必须经过一段时间的测试.第一个违规实例定义了带宽.我需要最大化这个带宽"/间隔.对于输出，我还需要知道哪个 L 导致这样的优化，因此我知道为给定的约束选择正确的 a.这就是正式的问题陈述.(希望这次我做对了)

Anyway, I now need to further optimize the problem with some constraints. To start, L is chosen. This allows me to calculate a. Once a is know, the polynomial is a known function of x only i.e p(x). I need to then determine the largest INTERVAL from 0->x over which the following constraint is satisfied: |dp(x)/dx - 1| < tol. This gives me a measure of the performance of the polynomial with the coefficient a. The interval is what I call the "bandwidth". I would like to emphasis two things: 1) The "bandwidth" is NOT the same as L. 2) All values of x within the "bandwidth" must meet the constraint. The function dp(x)/dx does oscillate in and out of the tolerance criteria, so testing the criteria for a single value of x does not work. It must be tested over an interval. The first instance of violation defines the bandwidth. I need to maximize this "bandwidth"/interval. For output, I also need to know which L lead to such an optimization, hence I know the correct a to choose for the given constraints. That is the formal problem statement. (I hope I got it right this time)

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

Now my problem is setting this whole thing up with MATLAB's optimization tools. I tried to follow ideas from the following articles:

• Optimization Toolbox™ 教程

为if 语句 设置optimize=1 将适用于约束优化.我想了一些嵌套优化是如何涉及的，但我无法进行任何工作.我从 IMSL 优化库中提供了该问题的已知解决方案以进行比较/检查.它们写在优化程序下面.无论如何，这是我到目前为止整理的代码:

Setting optimize=1 for the if statement will work with the constrained optimization. I thought some how nested optimization is involved, but I couldn't get anything to work. I provided known solutions to the problem from the IMSL optimization library to compare/check with. They are written below the optimization routine. Anyway, here is the code I've put together so far:

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 这样做过.我还应该注意，我上面的内容不起作用，并且对于受约束的优化是不正确的.

I could really use some help setting up the constrained optimization. I'm not only new to optimizations, I've never used MATLAB to do so. I should also note that what I have above does not work and is incorrect for the constrained optimization.

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

UPDATE: I added a for loop in the section if optimizeto show what I'm trying to achieve with the optimization. Obviously, I could just use this, but it seems very inefficient, especially if I increase the resolution of range and have to run this optimization many times. If you uncomment the plots, it will show how the bandwidth behaves. By looping over the full range, I'm basically testing every L but surely there's got to be a more efficient way to do this??

更新:已解决

推荐答案

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

So it seems fmincon is not the only tool for this job. In fact I couldn't even get it to work. Below, fmincon gets "stuck" on the IC and refuses to do anything...why...that's for a different post! Using the same layout and formulation, fminbnd finds the correct solution. The only difference, as far as I know, is that the former was using a conditional. But my conditional is nothing fancy, and really unneeded. So it's got to have something to do with the algorithm. I guess that's what you get when using a "black box". Anyway, after a long, drawn out, painful, learning experience, here is a solution:

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 的所有支持.如果没有帮助，我是不会想出来的！

I would like to especially thank @AndrasDeak for all their support. I wouldn't have figured it out without the assistance!

这篇关于如何使用匿名函数优化 MATLAB 中的约束积分表达式?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！