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

2023年2月5日 263次阅读

我有一个集成的错误表达式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,给出三个系数.

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

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

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

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

> Tutorial for the Optimization Toolbox™

为if语句设置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




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

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

更新：解决了

最佳答案 因此,似乎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的所有支持.没有帮助,我不会想出来的！

点赞