MATLAB Answers

AdarG
0

Find constraints on polynomial coefficients optimization

Asked by AdarG
on 11 Jul 2019
Latest activity Commented on by AdarG
on 13 Jul 2019
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?

  6 Comments

I have constraints on a1 and a2 in the form of
lb<a1*T^2+a2*T+a3<ub;
Those are not real constraints on the variables, only on theta.
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.

Sign in to comment.

3 Answers

Answer by Bruno Luong
on 12 Jul 2019
Edited by Bruno Luong
on 12 Jul 2019
 Accepted Answer

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

  5 Comments

% 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
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.
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?

Sign in to comment.


Answer by Matt J
on 11 Jul 2019
Edited by Matt J
on 11 Jul 2019

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.

  6 Comments

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.
My cost function throws an error when input out-of-bounds values, so I am forced to do the optimization in one go. :(
You should just set out of bound values to Inf.

Sign in to comment.


Answer by Matt J
on 11 Jul 2019
Edited by Matt J
on 11 Jul 2019

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

  0 Comments

Sign in to comment.