sqp() with trapz()

classic Classic list List threaded Threaded
17 messages Options
Reply | Threaded
Open this post in threaded view
|

sqp() with trapz()

dg.aragones
Hello all,

I have been working for a while in a code in Octave. All lines seem to work,
except for the last one, when I run an optimization using sqp(). This is the
code:

%Use Gnuplot for graphics in Octave
graphics_toolkit("gnuplot")

%Exact solution
F=@(mu) mu*log((-1-sqrt(1+mu^2))/(1-sqrt(1+mu^2)))-1; %Function for the
solver
mu=fzero(F,[0.01,1]) %Solution via Brent-Dekker-van Wijngaarden method
s=linspace(0,2,N=51); %s discretization
y=-sqrt(1+mu^2)+sqrt((s-1).^2+mu^2); %y exact values
%plot(s,y) %Exact solution

%Approximate solution
n=50; %Number of subintervals
h=2/n; %Step size
F=@(y) trapz(s,y); %Approximate objective function (Trapezoidal rule)
D=zeros(n+1,n+1);
D1=-eye(n/2)+diag(ones([n/2-1,1]),1);
D4=eye(n/2)+diag(-ones([n/2-1,1]),-1);
D(1:n/2,1:n/2)=D1;
D(end-n/2+1:end,end-n/2+1:end)=D4;
D(n/2,n/2+1)=1;
D(n/2+1,n/2+2)=1/2;
D(n/2+2,n/2+1)=-1;
D(n/2+1,n/2)=-1/2;
D=1/h*D; %Differentiation matrix
yp=@(y) transpose(D*transpose(y)); %Derivative (Finite Differences)
G=@(y) sqrt(1-yp(y).^2); %Function in the integral constraint
R=@(y) trapz(s,G(y))-1; %Approximation of the integral constraint
(Trapezoidal rule)
sqp(y,F,R) %Optimization via sequential quadratic programming

And the error it returns:

error: trapz: X and Y must have same shape
error: called from
    trapz at line 120 column 9
    sqp at line 351 column 7
    Catenary at line 28 column 1

However the trapz() function alone seems to work correctly. Does anyone know
what is happening here?

Thank you in advance for any help!



--
Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html


Reply | Threaded
Open this post in threaded view
|

Re: sqp() with trapz()

karl
Am 25.08.19 um 17:18 schrieb dg.aragones:

> Hello all,
>
> I have been working for a while in a code in Octave. All lines seem to work,
> except for the last one, when I run an optimization using sqp(). This is the
> code:
>
> %Use Gnuplot for graphics in Octave
> graphics_toolkit("gnuplot")
>
> %Exact solution
> F=@(mu) mu*log((-1-sqrt(1+mu^2))/(1-sqrt(1+mu^2)))-1; %Function for the
> solver
> mu=fzero(F,[0.01,1]) %Solution via Brent-Dekker-van Wijngaarden method
> s=linspace(0,2,N=51); %s discretization
> y=-sqrt(1+mu^2)+sqrt((s-1).^2+mu^2); %y exact values
> %plot(s,y) %Exact solution
>
> %Approximate solution
> n=50; %Number of subintervals
> h=2/n; %Step size
> F=@(y) trapz(s,y); %Approximate objective function (Trapezoidal rule)
> D=zeros(n+1,n+1);
> D1=-eye(n/2)+diag(ones([n/2-1,1]),1);
> D4=eye(n/2)+diag(-ones([n/2-1,1]),-1);
> D(1:n/2,1:n/2)=D1;
> D(end-n/2+1:end,end-n/2+1:end)=D4;
> D(n/2,n/2+1)=1;
> D(n/2+1,n/2+2)=1/2;
> D(n/2+2,n/2+1)=-1;
> D(n/2+1,n/2)=-1/2;
> D=1/h*D; %Differentiation matrix
> yp=@(y) transpose(D*transpose(y)); %Derivative (Finite Differences)
> G=@(y) sqrt(1-yp(y).^2); %Function in the integral constraint
> R=@(y) trapz(s,G(y))-1; %Approximation of the integral constraint
> (Trapezoidal rule)
> sqp(y,F,R) %Optimization via sequential quadratic programming
>
> And the error it returns:
>
> error: trapz: X and Y must have same shape
> error: called from
>      trapz at line 120 column 9
>      sqp at line 351 column 7
>      Catenary at line 28 column 1
>
> However the trapz() function alone seems to work correctly. Does anyone know
> what is happening here?
>
> Thank you in advance for any help!
>
>
>
> --
> Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html
>
>
Did you look up the description of trapz?

https://octave.sourceforge.io/octave/function/trapz.html




Reply | Threaded
Open this post in threaded view
|

Re: sqp() with trapz()

Doug Stewart-4
In reply to this post by dg.aragones


On Sun, Aug 25, 2019 at 11:18 AM dg.aragones <[hidden email]> wrote:
Hello all,

I have been working for a while in a code in Octave. All lines seem to work,
except for the last one, when I run an optimization using sqp(). This is the
code:

%Use Gnuplot for graphics in Octave
graphics_toolkit("gnuplot")

%Exact solution
F=@(mu) mu*log((-1-sqrt(1+mu^2))/(1-sqrt(1+mu^2)))-1; %Function for the
solver
mu=fzero(F,[0.01,1]) %Solution via Brent-Dekker-van Wijngaarden method
s=linspace(0,2,N=51); %s discretization
y=-sqrt(1+mu^2)+sqrt((s-1).^2+mu^2); %y exact values
%plot(s,y) %Exact solution

%Approximate solution
n=50; %Number of subintervals
h=2/n; %Step size
F=@(y) trapz(s,y); %Approximate objective function (Trapezoidal rule)
D=zeros(n+1,n+1);
D1=-eye(n/2)+diag(ones([n/2-1,1]),1);
D4=eye(n/2)+diag(-ones([n/2-1,1]),-1);
D(1:n/2,1:n/2)=D1;
D(end-n/2+1:end,end-n/2+1:end)=D4;
D(n/2,n/2+1)=1;
D(n/2+1,n/2+2)=1/2;
D(n/2+2,n/2+1)=-1;
D(n/2+1,n/2)=-1/2;
D=1/h*D; %Differentiation matrix
yp=@(y) transpose(D*transpose(y)); %Derivative (Finite Differences)
G=@(y) sqrt(1-yp(y).^2); %Function in the integral constraint
R=@(y) trapz(s,G(y))-1; %Approximation of the integral constraint
(Trapezoidal rule)
sqp(y,F,R) %Optimization via sequential quadratic programming

And the error it returns:

error: trapz: X and Y must have same shape
error: called from
    trapz at line 120 column 9
    sqp at line 351 column 7
    Catenary at line 28 column 1

However the trapz() function alone seems to work correctly. Does anyone know
what is happening here?

Thank you in advance for any help!



--
Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html



You have 2 different functions named F
which one do you want in the aqp function?


--
DASCertificate for 206392



Reply | Threaded
Open this post in threaded view
|

Re: sqp() with trapz()

dg.aragones
In reply to this post by karl
Thank you for your answer!

Yes, I did. However, I think I've written the code accordingly...



--
Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html


Reply | Threaded
Open this post in threaded view
|

Re: sqp() with trapz()

dg.aragones
In reply to this post by Doug Stewart-4
Thank you for your answer!

It's true. I want to use the second one in the sqp() function.



--
Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html


Reply | Threaded
Open this post in threaded view
|

Re: sqp() with trapz()

karl
In reply to this post by dg.aragones
Am 26.08.19 um 01:06 schrieb dg.aragones:

> Thank you for your answer!
>
> Yes, I did. However, I think I've written the code accordingly...
>
>
>
> --
> Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html
>
>
The error message says that

        error: trapz: X and Y must have same shape

You have written:

R=@(y) trapz(s,G(y))-1

Do s and G(y) have the same shape?



Reply | Threaded
Open this post in threaded view
|

Re: sqp() with trapz()

dg.aragones
Yes, they have the same size. It's weird...



--
Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html


Reply | Threaded
Open this post in threaded view
|

Re: sqp() with trapz()

karl
Am 26.08.19 um 18:25 schrieb dg.aragones:
> Yes, they have the same size. It's weird...
>
>
>
> --
> Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html
>
>
Please look at your code. Trapz expects two vectors of the same size.

You have the placeholder variable s (replace it by x since it can any name)

and then a function call. Are these two vectors of the same size?



Reply | Threaded
Open this post in threaded view
|

Re: sqp() with trapz()

karl
Am 26.08.19 um 19:07 schrieb [hidden email]:

> Am 26.08.19 um 18:25 schrieb dg.aragones:
>> Yes, they have the same size. It's weird...
>>
>>
>>
>> --
>> Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html
>>
>>
> Please look at your code. Trapz expects two vectors of the same size.
>
> You have the placeholder variable s (replace it by x since it can any name)
>
> and then a function call. Are these two vectors of the same size?
>
>
>
Replace s by y in the answer.



Reply | Threaded
Open this post in threaded view
|

Re: sqp() with trapz()

dg.aragones
Yes, they have exactly the same size:

s is a vector 51x1
G(y) returns a vector 51x1



--
Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html


Reply | Threaded
Open this post in threaded view
|

Re: sqp() with trapz()

karl
Am 26.08.19 um 20:12 schrieb dg.aragones:

> Yes, they have exactly the same size:
>
> s is a vector 51x1
> G(y) returns a vector 51x1
>
>
>
> --
> Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html
>
>
if you mean, I give up.



Reply | Threaded
Open this post in threaded view
|

Re: sqp() with trapz()

Octave - General mailing list
In reply to this post by dg.aragones
On 8/26/19 2:12 PM, dg.aragones wrote:
> Yes, they have exactly the same size:
>
> s is a vector 51x1
> G(y) returns a vector 51x1
>
The error happens in the other invocation of trapz, via the F function

F=@(y) trapz(s,y);



Reply | Threaded
Open this post in threaded view
|

Re: sqp() with trapz()

dg.aragones
Yes, but I don't understand why, because they have also the same size. They
are both 51x1 vectors.



--
Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html


Reply | Threaded
Open this post in threaded view
|

Re: sqp() with trapz()

Octave - General mailing list
On 8/26/19 2:32 PM, dg.aragones wrote:
> Yes, but I don't understand why, because they have also the same size. They
> are both 51x1 vectors.

Nope, one is 1x51 and the other is 51x1.

By using @ functions, it's harder to debug this because you can't put
tests/printouts into them, so I redefined them to use the function
trapz_debug, and defined it like so:

function f = trapz_debug (a, b)
   size (a)
   size (b)
   f = trapz (a, b)
endfunction




Reply | Threaded
Open this post in threaded view
|

Re: sqp() with trapz()

dg.aragones
Ok, thank you very much for your answer!



--
Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html


Reply | Threaded
Open this post in threaded view
|

Re: sqp() with trapz()

nrjank
In reply to this post by Octave - General mailing list
On Mon, Aug 26, 2019 at 2:39 PM Przemek Klosowski via Help list for GNU Octave <[hidden email]> wrote:
On 8/26/19 2:32 PM, dg.aragones wrote:

I copied and pasted your original code.  looking at the sqp line, I see in the help that sqp of the form:
sqp (X0, PHI, G)
requires:
X0 to be a vector
PHI to be the objective function that accepts a vector and returns a scalar and
G to be a constraint function that must accept a vector argument and return a vector.

If your case,  both F and R accept a vector and returns a scalar.  could that be part of the problem?  it gives the error: 

>> sqp(y,F,R)
error: operator *: nonconformant arguments (op1 is 51x51, op2 is 1x51)

(both s and y are 1x51 when I run your code unmodified.)

is sqp performing some array functions internally that is expecting a vector output from R to build a 51x51 array, but since it's only outputting scalars op2 is the wrong size?


Reply | Threaded
Open this post in threaded view
|

Re: sqp() with trapz()

dg.aragones
I have tried modifying the code to make R(y) return a vector, but it gives
the same error...



--
Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html