Find constraints on polynomial coefficients optimization

I am trying to find the optimal coefficients of the polynomial of the form:

theta=a1*t^2 +a2*t+a3 (i.e., to find a1,a2,a3) for some cost function.

I'm using patternsearch and I need to formulate the nonlinear/linear constraints on a1,a2,a3.

The problem is that I have constraints on theta (say [lb,ub]) and the range of t (say [0,T]), but not on the coefficients themselves.

So far, I've managed to formulate these constraints:

lb<a3<ub;
lb<a1*T^2+a2*T+a3<ub;

I can't figure out the 3rd constraint on the extrimum on t=-a2/(2*a1). I care only if is in the rectancle [0,T],[lb,ub].

Any idea?

[Walter Roberson]

How would you know if a set of coefficients were optimal?

[Matt J]

AdarG's comment relocated here:

Hi Walter,

patternsearch creates the a1,a2,a3 coefficients and use the polynomial to obtain a value of the cost function and minimize it. I only need to construct the constraints on the coefficients for the optimizator.

[Walter Roberson]

Is it correct that you are trying to find a1 a2 a3 t values that minimize theta under constraints on theta, t, and a3, but no constraints on a1 or a2?

If so then you can always choose theta as the minimum permitted: you have enough freedom to choose a1 and a2 to give any output you want.

[AdarG]

I have constraints on a1 and a2 in the form of

lb<a1*T^2+a2*T+a3<ub;

[Walter Roberson]

Those are not real constraints on the variables, only on theta.

[AdarG]

Now I understand. I will clarify,

1. I don't want to find the minimum of theta, but a minimum of some cost function that the polynomial is producing (it envolves a very complicated ode).
2. I have constraints on theta.
3. Since my design parameters are a1,a2,a3, I need to reformulate the constraints on them.

Answer 1

[Bruno Luong]

Why can't you implement

ts := max(min(-a2/(2*a1),T),0);

Then add the 6 constraints into your minimization pb:

two non-linear (and not differentiable):

lb <= theta(ts) <= ub

four non equality linear contstraints;

lb <= theta(0) <= ub
lb <= theta(T) <= ub

[Torsten]

... or just in nonlcon:

c = [];
ceq = [];
ts = -a2/(2*a1);
if  ts > 0 & ts < T
  c(1) = theta(ts) - ub;
  c(2) = lb - theta(ts);
end

plus the four inequality constraints that can be specified in A and b.

[AdarG]

Ok, I tried this and I think I am almost there.

I get an error:

Unable to perform assignment because the left and right sides have a different number of elements.

Any idea why this error appears?

Some details:

The error originates in line 155 in augLagFun.m

[cin(:),ceq(:)] = feval(nonlcon,points,conFcnArg{:});

My nonlinconstraint function is defined as a nested function within the same function that calls the optimizator and only when c(1) and c(2) have a value, the error appears.

apart from design_param, all other params are global within the nested functions.

design_param is actually a vecotor of 5 elements; 3 of them are used for the polynomial of theta; the other 2 are constant values of other variables to be optimized.

function [c,ceq] = theta_extremum_con(design_param)
        c=[];
        if design_param(1)~=0 %if a1 is zero, theta is linear so no extremum
            ts = -design_param(2)/(2*design_param(1));
            if  ts > indep_var_bounds(1) && ts < indep_var_bounds(2)
                c(1) = polyval(design_param(1:3),ts) - extrm_val(2);
                c(2) = extrm_val(1) - polyval(design_param(1:3),ts);
            end
        end
        ceq = [];
end

[Bruno Luong]

% Generate test data
T = 1;
a = randn(3,1);
tin = linspace(0,T,20);
thetafun = @(a,t) polyval(a(1:3),t);
yin = thetafun(a,tin);
yin = yin + 0.03*randn(size(yin));
lb = -1;
ub = 1;

% Fit parabola to (ti,yi) with constraints
a0 = zeros(3,1);
A = [0; T].^(2:-1:0);
A = [A; -A];
b = repelem([ub; -lb],2);
fun = @(a) l2(a,tin,yin,thetafun);
nonlcon = @(a) nlcon(a,T,lb,ub,thetafun);
a = fmincon(fun,a0,A,b,[],[],[],[],nonlcon);

% Check with graphics
close all
h1 = plot(tin,yin,'.r');
hold on
ti = linspace(0,T,100);
yi = thetafun(a,ti);
h3 = plot([0 0; T T], [lb ub; lb ub],'g');
h2 = plot(ti,yi,'b');
legend([h1,h2,h3(1)],'data','fit','limits');

%%
function f = l2(a,tin,yin,thetafun)
f = sum((thetafun(a,tin)-yin).^2);
end

%%
function [c,ceq] = nlcon(a,T,lb,ub,thetafun)
ts = max(min(-a(2)/(2*a(1)),T),0);
yts = thetafun(a,ts);
c = [yts-ub;
     lb-yts];
ceq = [];
end

[AdarG]

Thanks, that did the trick.

I used your function

function [c,ceq] = nlcon(a,T,lb,ub,thetafun)
ts = max(min(-a(2)/(2*a(1)),T),0);
yts = thetafun(a,ts);
c = [yts-ub;
     lb-yts];
ceq = [];
end

I understood that c must give the same number of elements at every call, so I can't use the syntax I was using with the if statement.

[AdarG]

Although this issue is solved, when I run the optimization, the optimizator runs alot of function evaluations but with no iteration progress, so the optimization is verrrrrry slow.

I checked that the non linear constratint function works properly and it does.

Any idea why so many evaluations performed?

Answer 2

[Matt J]

What's to figure out? You've already articulated that the (nonlinear) constraints on the extremum are,

0<=-a2/(2*a1)<=T

The only thing I might recommend is that converting them to linear constraints,

0<=-a2<=2*T*a1
a1>=0

might make things easier for patternsearch.

[AdarG]

Hi Matt,

I don't mind the extremum location to be in the range

0<=-a2/(2*a1)<=T

I mind that if it is in the range, then the extremum must not exceed the boundaries on theta.

[Matt J]

I see. Well, then you might divide the optimization into 3 sub-problems, corresponding to the 3 different regions where the critical t can lie. For each sub-problem, you apply a different set of constraints:

Case I constraints:

lb<=a3<=ub;
lb<= a1*T^2+a2*T+a3 <=ub;
-a2/(2*a1)<=0

Case II constraints:

lb<=a3<=ub;
lb<= a1*T^2+a2*T+a3 <=ub;
-a2/(2*a1)>=T

Case III constraints:

lb<=a3<=ub;
lb<= a1*T^2+a2*T+a3 <=ub;
0<=-a2/(2*a1)<=T
lb<=polyval(a,-a2/(2*a1))<=ub

[AdarG]

Thanks Matt, but how do I implement it practically? Should I call the optimization routine 3 times?

If I call only one time, how do I construct the A matrix and b vector?

In the meantime, I implement the upper 2 constraints ( lb<=a3<=ub;

lb<= a1*T^2+a2*T+a3 <=ub;) and check in the optimization (after the design parameters were picked) if the design parameters meet the extremum condition.

[Matt J]

Thanks Matt, but how do I implement it practically? Should I call the optimization routine 3 times?

Yes, you solve each sub-case and take the case with the best optimal value.

If I call only one time, how do I construct the A matrix and b vector?

Each row of the matrix corresponds to one of the linear constraints.

[AdarG]

My cost function throws an error when input out-of-bounds values, so I am forced to do the optimization in one go. :(

[Matt J]

You should just set out of bound values to Inf.

Answer 3

[Matt J]

It might be more natural here to use fseminf,

x = fseminf(fun,[a1,a2,a3], 2,  @(a,s) seminfcon(a,s,T,lb,ub));


function [c,ceq,K_ub,K_lb,s]= seminfcon(a,s,T,lb,ub)

    % No finite nonlinear inequality and equality constraints
    c = [];
    ceq = [];
    
    % Sample set
    if isnan(s(1))
        % Initial sampling interval
        s = [0.01 0; 0.01 0];
    end
    
    t1 = 0:s(1):T;
    t2 = 0:s(2):T;

    
    % Evaluate the semi-infinite constraint
    K_ub = polyval(a,t1)-ub;
    K_lb = lb - polyval(a,t2);
end