Thread Subject:

Find min/max analiticaly

Hello,

I'm working in an optimization problem with SNOPT. I'm trying to change the control discretization from a linear one to higher order. I implemented the PCHIP but I would like to implement also the SPLINE but the problem with the bplines is that we have overshoot between the points and because of that I'll have to change all the Constrains parameters (because now the max and min values could be bigger/lower than a point).
The problem is I can't use Max(A) because I will need the derivative of the constrains functions later so I have to find a way to obtain the min and max values analytically.

Thank you

"Marc " <job21marc@hotmail.com> wrote in message <jplbcp$9ud$1@newscl01ah.mathworks.com>...
>
> The problem is I can't use Max(A) because I will need the derivative of the constrains functions later so I have to find a way to obtain the min and max values analytically.
=======

It's highly doubtful you'll be able to do this.

"Matt J" wrote in message <jplrkb$pr0$1@newscl01ah.mathworks.com>...
> "Marc " <job21marc@hotmail.com> wrote in message <jplbcp$9ud$1@newscl01ah.mathworks.com>...

>
> It's highly doubtful you'll be able to do this.

In theory it's possible.

A cubic spline has the second derivative at the knots computed as solution of a triangular matrix. The triangular matrix depend only on knot location. So the second derivative is of the form

[f"(x_i)] = A * y_i

where A is (n x n) matrix and y_i is f(x_i).

The local optima of f are where f(x) = 0, they can be then computed by solving a piecewise second order polynomial, which can be expressed analytic as f"(x_i) and f(x_i).

Then just evaluate f at local optima, and take a max of those.

There might be a bound formula of max than can be even more easily computed from the 2-norm of A.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <jpm3ug$5d4$1@newscl01ah.mathworks.com>...
> "Matt J" wrote in message <jplrkb$pr0$1@newscl01ah.mathworks.com>...
> > "Marc " <job21marc@hotmail.com> wrote in message <jplbcp$9ud$1@newscl01ah.mathworks.com>...
>
> >
> > It's highly doubtful you'll be able to do this.
>
> In theory it's possible.
>
> A cubic spline has the second derivative at the knots computed as solution of a triangular matrix. The triangular matrix depend only on knot location. So the second derivative is of the form
>
> [f"(x_i)] = A * y_i
>
> where A is (n x n) matrix and y_i is f(x_i).
>
> The local optima of f are where f(x) = 0, they can be then computed by solving a piecewise second order polynomial, which can be expressed analytic as f"(x_i) and f(x_i).
>
> Then just evaluate f at local optima, and take a max of those.
>
> There might be a bound formula of max than can be even more easily computed from the 2-norm of A.
>
> Bruno

wow, thank you Bruno but I don't see how to make it. my code is simple,
for example, how can you calculated min and max from that:
>> x=1:10;
>> x=0:10;
>> y=[3 4 8 9 5 7 2 9 4 6 1];
>> cs=spline(x,y);

Thanks

"Marc " <job21marc@hotmail.com> wrote in message <jpngc1$r4c$1@newscl01ah.mathworks.com>...

>
> wow, thank you Bruno but I don't see how to make it. my code is simple,
> for example, how can you calculated min and max from that:
> >> x=1:10;
> >> x=0:10;
> >> y=[3 4 8 9 5 7 2 9 4 6 1];
> >> cs=spline(x,y);
>

x=0:10;
y=[3 4 8 9 5 7 2 9 4 6 1];
cs=spline(x,y);

coefs = cs.coefs;
a = 3*coefs(:,1);
b = 2*coefs(:,2);
c = 1*coefs(:,3);
delta = b.^2 - 4*a.*c;
x1 = (-b - sqrt(delta)) ./ (2*a);
x2 = (-b + sqrt(delta)) ./ (2*a);
left = cs.breaks';
h = diff(left);
b1 = delta >= 0 & x1 >= 0 & x1 <= h;
b2 = delta >= 0 & x2 >= 0 & x2 <= h;
x = [left(1);
       left(b1) + x1(b1); ...
       left(b2) + x2(b2);
      left(end)];
 
[ymax imax] = max(ppval(cs,x));
xmax = x(imax);

fprintf('ymax = %g\n', ymax);
xi = linspace(0,10);
figure;
plot(xi, ppval(cs, xi), 'b', xmax, ymax, 'or')

% Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <jpo84r$8ao$1@newscl01ah.mathworks.com>...
> "Marc " <job21marc@hotmail.com> wrote in message <jpngc1$r4c$1@newscl01ah.mathworks.com>...
>
> >
> > wow, thank you Bruno but I don't see how to make it. my code is simple,
> > for example, how can you calculated min and max from that:
> > >> x=1:10;
> > >> x=0:10;
> > >> y=[3 4 8 9 5 7 2 9 4 6 1];
> > >> cs=spline(x,y);
> >
>
> x=0:10;
> y=[3 4 8 9 5 7 2 9 4 6 1];
> cs=spline(x,y);
>
> coefs = cs.coefs;
> a = 3*coefs(:,1);
> b = 2*coefs(:,2);
> c = 1*coefs(:,3);
> delta = b.^2 - 4*a.*c;
> x1 = (-b - sqrt(delta)) ./ (2*a);
> x2 = (-b + sqrt(delta)) ./ (2*a);
> left = cs.breaks';
> h = diff(left);
> b1 = delta >= 0 & x1 >= 0 & x1 <= h;
> b2 = delta >= 0 & x2 >= 0 & x2 <= h;
> x = [left(1);
> left(b1) + x1(b1); ...
> left(b2) + x2(b2);
> left(end)];
>
> [ymax imax] = max(ppval(cs,x));
> xmax = x(imax);
>
> fprintf('ymax = %g\n', ymax);
> xi = linspace(0,10);
> figure;
> plot(xi, ppval(cs, xi), 'b', xmax, ymax, 'or')
>
> % Bruno

Thanks so much Bruno, now i stil have a problem that i don't know how to solve....
I'm doing two interpolation, this code works perfect with one but with the other is more complicated to implement i think.

The point is that I interpolate a spline([157x1],[168x157]) and the result is:

      form: 'pp'
    breaks: [1x157 double]
     coefs: [26208x4 double]
    pieces: 156
     order: 4
       dim: 168
I saw that they put all the coefs together but this is not the problem because I can do it with a for and step-by-step. The problem is that now x1 and x2 are 168x1 and h is 156X1 and I have a problem in:
b1 = delta >= 0 & x1 >= 0 & x1 <= h;
b2 = delta >= 0 & x2 >= 0 & x2 <= h;
because they are not of the same size....

the 168X157 is a diagonal matrix of 1 in the diagonal and 0 rest and the 157x1 is a 0:0.02:3 and adding all the 0.25-0.75-1.25-1.75-2.25-2.75.

thanks for everything again, you're helping me so much.

Multiple 1d spline curves:

x=0:10;
y=rand(5,length(x));
cs=spline(x,y);

xi = linspace(0,10);
figure;
plot(xi, ppval(cs, xi));
hold on;

n = size(y,1);
left = cs.breaks';
h = diff(left);
for k=1:size(y,1)
    csk = cs;
    csk.dim = 1;
    coefs = cs.coefs(k:n:end,:);
    csk.coefs = coefs;
    a = 3*coefs(:,1);
    b = 2*coefs(:,2);
    c = 1*coefs(:,3);
    delta = b.^2 - 4*a.*c;
    sqrt_d = sqrt(delta);
    x1 = (-b - sqrt_d) ./ (2*a);
    x2 = (-b + sqrt_d) ./ (2*a);
    
    b1 = delta >= 0 & x1 >= 0 & x1 <= h;
    b2 = delta >= 0 & x2 >= 0 & x2 <= h;
    xmax = [left(1);
            left(b1) + x1(b1); ...
            left(b2) + x2(b2);
            left(end)];
    
    [ymax imax] = max(ppval(csk,xmax));
    xmax = xmax(imax);
    
    fprintf('ymax = %g\n', ymax);
    plot(xmax, ymax, 'or');
end

% Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <jpqklr$ekv$1@newscl01ah.mathworks.com>...
> Multiple 1d spline curves:
>
> x=0:10;
> y=rand(5,length(x));
> cs=spline(x,y);
>
> xi = linspace(0,10);
> figure;
> plot(xi, ppval(cs, xi));
> hold on;
>
> n = size(y,1);
> left = cs.breaks';
> h = diff(left);
> for k=1:size(y,1)
> csk = cs;
> csk.dim = 1;
> coefs = cs.coefs(k:n:end,:);
> csk.coefs = coefs;
> a = 3*coefs(:,1);
> b = 2*coefs(:,2);
> c = 1*coefs(:,3);
> delta = b.^2 - 4*a.*c;
> sqrt_d = sqrt(delta);
> x1 = (-b - sqrt_d) ./ (2*a);
> x2 = (-b + sqrt_d) ./ (2*a);
>
> b1 = delta >= 0 & x1 >= 0 & x1 <= h;
> b2 = delta >= 0 & x2 >= 0 & x2 <= h;
> xmax = [left(1);
> left(b1) + x1(b1); ...
> left(b2) + x2(b2);
> left(end)];
>
> [ymax imax] = max(ppval(csk,xmax));
> xmax = xmax(imax);
>
> fprintf('ymax = %g\n', ymax);
> plot(xmax, ymax, 'or');
> end
>
> % Bruno

Thank you so much Bruno! you helped me a lot! I'm gonna try to implement it, I find anything else I will post it

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message
> xmax = [left(1);
> left(b1) + x1(b1); ...
> left(b2) + x2(b2);
> left(end)];

> % Bruno

Bruno thank you in advance again but another question,

Do you think that we can be able to differentiate the xmax equation?

Because I need the differential equation of xmax to be able to find the max and min when my values changes.

thank you so much again

Marc

"Marc " <job21marc@hotmail.com> wrote in message <jq7crh$pqt$1@newscl01ah.mathworks.com>...

> Do you think that we can be able to differentiate the xmax equation?
>

Yes. Not very straightforward to derive, but possible.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <jq7f4d$4ua$1@newscl01ah.mathworks.com>...
> "Marc " <job21marc@hotmail.com> wrote in message <jq7crh$pqt$1@newscl01ah.mathworks.com>...
>
> > Do you think that we can be able to differentiate the xmax equation?
> >
>
> Yes. Not very straightforward to derive, but possible.
>
> Bruno

It will be great if you could help me with that or tell me some help.

Thank you

"Marc " <job21marc@hotmail.com> wrote in message <jqkmq6$9sr$1@newscl01ah.mathworks.com>...

>
> It will be great if you could help me with that or tell me some help.

This task is too big, and I can't afford for that.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <jqku0m$9ae$1@newscl01ah.mathworks.com>...
> "Marc " <job21marc@hotmail.com> wrote in message <jqkmq6$9sr$1@newscl01ah.mathworks.com>...
>
> >
> > It will be great if you could help me with that or tell me some help.
>
> This task is too big, and I can't afford for that.
>
> Bruno

ok, no worries. You helped me a lot,

thank you so much