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.
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.
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.
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)
Now my problem is setting this whole thing up with MATLAB's optimization tools. I tried to follow ideas from the following articles:
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
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: 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??
UPDATE: Solved
