Newbie question - sqp solver equality constraints function

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
14 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Newbie question - sqp solver equality constraints function

alexegin
Hello all!
I am trying to use sqp solver to find optimal vector X. All of the vector elements must meet the following constraints:
# 1. All of the elements must be more than -1.0 and less than 1.0

-1.0 < X(i) < 1.0

# 2.1 Sum of all elements those are less than zero must be equal to -1.0
# 2.2 Sum of all the elements those are more than zero must be equal to 1.0

dMinuses = 0.0;
dPluses   = 0.0;

for i = 1:rows( X )
  if ( X( i ) < 0.0 )
    dMinuses += X( i );
  else
    dPluses += X( i );
  endif;
endfor;

bResult = ( ( dMinuses + 1.0 ) == 0.0 ) && ( ( dPluses - 1.0 ) == 0.0 );
The first constraint -1.0 < X(i) < 1.0 can be easily passed to sqp solver by upper and lower bounds. But I do not know how to pass the second ones simultaneously:
SUM( X( i ) ) - 1.0 = 0.0, X( i ) >= 0.0
SUM( X( i ) ) + 1.0 = 0.0, X( i ) < 0.0
Thank you very much in advance.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newbie question - sqp solver equality constraints function

Juan Pablo Carbajal-2
On Wed, Jun 21, 2017 at 6:40 PM, alexegin <[hidden email]> wrote:

> Hello all!
> I am trying to use *sqp* solver to find optimal vector *X*. All of the
> vector elements /must/ meet the following constraints:
> # 1. All of the elements must be more than -1.0 and less than 1.0-1.0 < X(i)
> < 1.0# 2.1 Sum of all elements those are _less than zero_ must be equal to
> -1.0# 2.2 Sum of all the elements those are _more than zero_ must be equal
> to 1.0dMinuses = 0.0;dPluses   = 0.0;for i = 1:rows( X )  if ( X( i ) < 0.0
> )    dMinuses += X( i );  else    dPluses += X( i );  endif;endfor;bResult =
> ( ( dMinuses + 1.0 ) == 0.0 ) && ( ( dPluses - 1.0 ) == 0.0 );
> The first constraint -1.0 < X(i) < 1.0 can be easily passed to *sqp* solver
> by upper and lower bounds.But I do not know how to pass the second ones
> *simultaneously*:
> SUM( X( i ) ) - 1.0 = 0.0, X( i ) >= 0.0SUM( X( i ) ) + 1.0 = 0.0, X( i ) <
> 0.0
> Thank you very much in advance.
>
Hi,
As you said the range of X should be passed using UB, LB arguments.

The others you can pass using the equality constraints

function E = eq_cons (X)
 E = [sum(X(X<0)); sum(X(X>0))] - [-1; 1];
endfunction

Note that this function is problematic because of the jump introduced
when an element changes signs, but give it a try.

I am not sure but this alternative might be more suited for sqp (but
problably slower due to the copy of X)

function E = eq_cons (X)
 Y = X;
 Y(Y<0) *= -1;
 E = sum([X Y]).' - [0; 2];
endfunction

I get the same results for dim(x) = 10,..,100 with little difference
sin perfomrnace, but maybe your cost is more complex than my test:

g =@(X)  [sum(X(X<0)); sum(X(X>0))] - [-1; 1];
c =@(x) sumsq(x);
x0 = 2*rand(100,1)-1;
tic; [x obj info iter nf l] = sqp(x0,c,g,[],-1,1); toc

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newbie question - sqp solver equality constraints function

alexegin
This post was updated on .
Hi, Juan Pablo Carbajal-2!
Thank you very much for answering me. Your test is working perfectly:
clear all;
g =@( x ) [ sum( x( x < 0.0 ) ); sum( x( x >= 0.0 ) ) ] - [ -1.0; 1.0 ];
c =@( x ) sumsq( x ); 
x0 = 2.0 * rand( 10, 1 ) - 1.0; 
tic;
[ x obj info iter nf l ] = sqp( x0, c, g, [], -1.0, 1.0 );
toc;
But when I am trying to use constraint function in my script I get the following error:
error: qp: equality constraint matrix must be full row rank
Here is my script:
clear all;

global R = dlmread( "r.csv", "\t" ); # double [ 10x1 ]
global Q = dlmread( "q.csv", "\t" ); # double [ 10x10 ]
global Y = dlmread( "y.csv", "\t" ); # double [ 10x1 ]

global lb = dlmread( "lb.csv", "\t" ); # double [ 10x1 ]
global ub = dlmread( "ub.csv", "\t" ); # double [ 10x1 ]

g = @( x ) [ sum( x( x < 0.0 ) ); sum( x( x >= 0.0 ) ) ] - [ -1.0; 1.0 ];

function z = foo( m, n )  z = 0; for i = 1:rows( m ) z = z + m( i, 1 ) * n( i, 1 ); endfor; endfunction;

function y = phi( x ) global Y; y = foo( Y, x ) * -1; endfunction

tic;
[ x obj info iter nf l ] = sqp( R, @phi, g, [], lb, ub );
toc;
What am I doing wrong? Please do not tell me everything
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newbie question - sqp solver equality constraints function

Juan Pablo Carbajal-2
On Thu, Jun 22, 2017 at 9:04 AM, alexegin <[hidden email]> wrote:
clear all;
global R = dlmread( "r.csv", "\t" ); # double [ 10x1 ]
global Q = dlmread( "q.csv", "\t" );# double [ 10x10 ]
global Y = dlmread( "y.csv", "\t" ); # double [ 10x1 ]
global lb = dlmread( "lb.csv", "\t" ); # double [10x1 ]
global ub = dlmread( "ub.csv", "\t" ); # double [ 10x1 ]
g = @( x ) [ sum( x( x < 0.0 ) ); sum( x( x >= 0.0 ) ) ] - [ -1.0; 1.0 ];
function z = foo( m, n )
  z = 0;
  for i = 1:rows( m )
    z = z + m( i, 1 ) * n( i, 1 );
 endfor;
endfunction;
function y = phi( x )
   global Y;
   y = foo( Y, x ) * -1;
endfunction
tic;
[ x obj info iter nf l ] = sqp( R, @phi, g, [], -1.0, 1.0);
toc;
> What am I doing wrong? Please do not tell me /everything/

ok, lets simplify...get use to do this (also try to avoid that I get
your code in a single line, makes it hard to read).

1. foo is just sum(m .* n), or m.'*n . Note that you need to
transpose, to avoid doing it every time do Y = Y.';
2. Then your cost is phi =@(x) - Y*x (so you are searching maximum
projection of x on Y. No need of global variables. Avoid them wen
possible!)
3. Your script works here with fictitious data

Y = randn(10,1).';
phi =@(x) - Y*x;
g=@(x) [sum(x(x<0.0)); sum(x(x>=0.0))] - [ -1.0; 1.0 ];
x0 = 2*rand(10,1)-1;
[x obj info iter nf l] = sqp( x0, phi, g, [], -1.0, 1.0);

possible causes of error that you need to check:
1. Make sure R (the initial guess) complies with your UB and LB
2. Make sure Y has meaningful values.
4. Remove all globals and do clear all.

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newbie question - sqp solver equality constraints function

alexegin
Hi, "Juan Pablo Carbajal-2"<br/>

Juan Pablo Carbajal-2 wrote
... also try to avoid that I get your code in a single line, makes it hard to read
Sorry, I format my code using pre html tags - I am just trying to make it easier to read and it looks fine in a browser. But if it makes it hard to read I will not use HTML format.
I think I know what is the reason of the error - when I was checking intitial guess as you sad I realized that it contains positive values only. When at least one of the values is negative - everything runs without any errors. I do not know yet how to avoid it, but thank you very much !!!
Can I ask you one more time to help me with the Octave syntax  - as far as you can see I am not familiar with it
Here is my full script and I think moo function can be implemented in a much more easier way - simplified, as you said:

clear all;

function r = moo( A, B )
 
  Tmp = zeros( 1, rows( B ) ); # Temp matrix

  Bt = B'; # Transpose B
 
  # Fill in temp matrix
  for i = 1:columns( Bt )
    for j = 1:rows( A )
      Tmp( 1, i ) += A( j, i ) * Bt( 1, j );
    endfor
  endfor
 
  # Compute result
  #r = 0;
  #for i = 1:columns( Tmp )
  #  r = r + Tmp( 1, i ) * B( i, 1 );
  #endfor
 
  r = Tmp * B; # I think it is better, right?

endfunction

Y = dlmread( "y.csv", "\t" );
Yt = Y';

R = dlmread( "r.csv", "\t" );

Q = dlmread( "q.csv", "\t" );

f = dlmread( "f.csv", "\t" );

Lb = dlmread( "lb.csv", "\t" ); # -1.0
Ub = dlmread( "ub.csv", "\t" ); # 1.0
 
phi = @( x ) - Yt * x;
 
g = @( x ) [ sum( x( x < 0.0 ) ); sum( x( x >=0.0 ) ) ] - [ Lb; Ub ];

h = @( x ) f - moo( Q, x );

tic;
[x obj info iter nf l] = sqp( R, phi, g, h, Lb, Ub );
toc;

Thank you very much again - you helped me a lot!
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newbie question - sqp solver equality constraints function

Juan Pablo Carbajal-2
On Thu, Jun 22, 2017 at 12:41 PM, alexegin <[hidden email]> wrote:
> Juan Pablo Carbajal-2 wrote
>> ... also try to avoid that I get your code in a single line, makes it hard
>> to read
>
> Sorry, I format my code using *pre* html tags - I am just trying to make it
> easier to read and it looks fine in a browser. But if it makes it hard to
> read I will not use HTML format.
Yes, use plain text.

> I think I know what is the reason of the error - when I was checking
> intitial guess as you sad I realized that it contains *positive values
> only*. When at least one of the values is *negative* - everything runs
> without any errors. I do not know yet how to avoid it, but *thank you very
> much !!!*
The initial guess is something kind of arbitrary, just use a hand
crafted initial guess. the better the guess the faster you will get
your solution.

> Can I ask you one more time to help me with the Octave syntax  - as far as
> you can see I am not familiar with it
Note:
1. If you are using objects with real values (i.e. no complex values)
use .' for transpose to avoid calling conj().
2. Think in linear algebra not in loops (you probably come from a
language that is good at loops), try to use operators whenever you
can.

> Here is my full script and I think *moo* function can be implemented in a
> much more easier way - /simplified/, as you said:
>
> clear all;
>
> function r = moo( A, B )
>
>   Tmp = zeros( 1, rows( B ) ); # Temp matrix
Put declarations close to the point where the object is used.
>
>   Bt = B'; # Transpose B
Use .' (unless there are complex values)

>
>   # Fill in temp matrix
>   for i = 1:columns( Bt )
>     for j = 1:rows( A )
>       Tmp( 1, i ) += A( j, i ) * Bt( 1, j );
>     endfor
>   endfor
>
Tmp is just
Tmp = B.' * A;
No need of allocation

>   # Compute result
>   #r = 0;
>   #for i = 1:columns( Tmp )
>   #  r = r + Tmp( 1, i ) * B( i, 1 );
>   #endfor
>
>   r = Tmp * B; # I think it is better, right?
This is correct so your whole function is just
r =@(A,b) b.' * A * b
or even better
h = @(x) f - x.' * A * x

>
> endfunction
>
> Y = dlmread( "y.csv", "\t" );
> Yt = Y';
>
> R = dlmread( "r.csv", "\t" );
>
> Q = dlmread( "q.csv", "\t" );
>
> f = dlmread( "f.csv", "\t" );
>
> Lb = dlmread( "lb.csv", "\t" ); # -1.0
> Ub = dlmread( "ub.csv", "\t" ); # 1.0
>
> phi = @( x ) - Yt * x;
>
> g = @( x ) [ sum( x( x < 0.0 ) ); sum( x( x >=0.0 ) ) ] - [ Lb; Ub ];
>
> h = @( x ) f - moo( Q, x );
>
> tic;
> [x obj info iter nf l] = sqp( R, phi, g, h, Lb, Ub );
> toc;
>
> Thank you very much again - you helped me a lot!
>
Your problem seems quite simple, check using QP (maybe you can render
it in that form) instead of SQP. Anyways SQP should work.

>
>
> --
> View this message in context: http://octave.1599824.n4.nabble.com/Newbie-question-sqp-solver-equality-constraints-function-tp4683834p4683861.html
> Sent from the Octave - General mailing list archive at Nabble.com.
>
> _______________________________________________
> Help-octave mailing list
> [hidden email]
> https://lists.gnu.org/mailman/listinfo/help-octave

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newbie question - sqp solver equality constraints function

alexegin
Hi, Juan Pablo Carbajal-2

Juan Pablo Carbajal-2 wrote
... or even better h = @(x) f - x.' * A * x
I tried to implement non-equality constraint function that way and the result is incredible - it is a hundred times faster than the one with the loop! One hundred, man! It took me so much time to answer because it was very hard to believe and I was checking the results.
Thank you very, very much!!!

Juan Pablo Carbajal-2 wrote
Your problem seems quite simple, check using QP (maybe you can render
it in that form) instead of SQP. Anyways SQP should work.
Maybe I would, but now it seems to me that sqp is doing quite good - I mean the speed :)
Thank you very much!
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newbie question - sqp solver equality constraints function

Juan Pablo Carbajal-2
> Thank you very, very much!!!
I am glad to help!
If you ever had the chance, drop a donation to Octave
(https://www.gnu.org/software/octave/doc/interpreter/How-You-Can-Contribute-to-Octave.html
or directly to John W. Eaton
http://jweaton.org/?page_id=48), are the developers who give us this
excellent piece of software.

Good luck and come back soon!

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newbie question - sqp solver equality constraints function

alexegin
This post was updated on .
Juan Pablo Carbajal-2 wrote
Good luck and come back soon!
Hello again :)
I have the following code working:

clear all;
Y = 2.0 * rand( 5, 1 ) - 1.0;
f = 0.0025;
Lb = -1.0; # Lower bound
Ub = 1.0; # Upper bound

# Initial guess
R = ones( rows( Y ), 1 );
R( 2:2:end ) = R( 2:2:end ) * ( 2.0 / rows( Y ) ) * Ub;
R( 1:2:end ) = R( 1:2:end ) * ( 2.0 / rows( Y ) ) * Lb;

# Objective function
phi = @( x ) sumsq( x );

# Equality constraint function
g = @( x ) [ sum( x( x <= 0.0 ) ); sum( x( x >= 0.0 ) ) ] - [ Lb; Ub ];

# Non-equality constraint function
h = @( x ) f - x.' * Y;

# Solver
[ x, obj, info, iter, nf, l ] = sqp( R, phi, g, h, Lb, Ub );

But if I set lower bound to zero or above ( Lb = 0.0 or Lb = 0.5 ) I have the following error:

error: qp: equality constraint matrix must be full row rank

What am I doing wrong? I am stuck ...
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newbie question - sqp solver equality constraints function

Juan Pablo Carbajal-2
> # Equality constraint function
> g = @( x ) [ sum( x( x <= 0.0 ) ); sum( x( x >= 0.0 ) ) ] - [ Lb; Ub ];
> # Solver
> [ x, obj, info, iter, nf, l ] = sqp( R, phi, g, h, Lb, Ub );
>
> But if I set lower bound to zero ( *Lb = 0.0* ) I have the following error:
>
> *error: qp: equality constraint matrix must be full row rank*
>
> What am I doing wrong? I am stuck ...

My guess is that, one of the components of your Equality constraint is
empty and hence its derivative is zero, giving zero eigenvalue in the
matrix assembled for QP.
Nevertheless the function g is not meaningful when Lb or Ub are 0. You
are saying that the sum of non-positive (non-negative) components
should be zero, if the components are not zero, this cannot be
fulfilled (sum of nonzero numbers of equal sign cannot be zero), hence
the only solution is all components to be zero. Do you set the initial
guess with this property?

What are you trying to achieve?

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newbie question - sqp solver equality constraints function

alexegin
Juan Pablo Carbajal-2 wrote
Nevertheless the function g is not meaningful when Lb or Ub are 0  ...
Do you set the initial guess with this property?
What are you trying to achieve?
You are right - when Lb ( Ub ) is zero, equality constraint function becomes just:

g = @( x ) sum( x ) - Ub;

I am trying to achieve a kind of universal solution - I run Octave script from another process, prepearing data for it (see dlmread above). I would like to have an opportunity to set Lb to -1.0 or to 0.0 while runnig the same script, somewhat:

#  I am not sure the following syntax is possible
if ( Lb < 0.0 )
    g = @( x ) [ sum( x( x <= 0.0 ) ); sum( x( x >= 0.0 ) ) ] - [ Lb; Ub ];
else
    g = @( x ) sum( x ) - Ub;
endif

As for initial guess, I set it this way (assume that Lb = -1.0):

R = ones( rows( Y ), 1 );
R( 2:2:end ) = R( 2:2:end ) * ( 2.0 / rows( Y ) ) * Ub; # sum( R( R >= 0.0 ) ) = 1.0
R( 1:2:end ) = R( 1:2:end ) * ( 2.0 / rows( Y ) ) * Lb; # sum( R( R <= 0.0 ) ) = -1.0
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newbie question - sqp solver equality constraints function

Juan Pablo Carbajal-2
> if ( Lb < 0.0 )
>     g = @( x ) [ sum( x( x <= 0.0 ) ); sum( x( x >= 0.0 ) ) ] - [ Lb; Ub ];
> else
>     g = @( x ) sum( x ) - Ub;
> endif
>
Yeah, that's fine.
Make sure your initial guess is also compatible with Lb and Ub.

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newbie question - sqp solver equality constraints function

alexegin
Juan Pablo Carbajal-2 wrote
...
Yeah, that's fine.
Make sure your initial guess is also compatible with Lb and Ub.
It is amazing, but  - yes, it is fine! After answering you I put it that way and - it works! Now I have the following code with the initial guess setting:

R = ones( rows( Y ), 1 );
if ( Lb < 0.0 )
  R( 2:2:end ) = R( 2:2:end ) * ( ( Ub - Lb ) / rows( Y ) ) * Ub;
  R( 1:2:end ) = R( 1:2:end ) * ( ( Ub - Lb ) / rows( Y ) ) * Lb;
  R( 2, 1 ) += 1.0 - sum( R( R >= 0.0 ) );
  R( 1, 1 ) -= sum( R( R <= 0.0 ) ) + 1.0;
else
  R = R * ( Ub + Lb ) / rows( Y );
endif

if ( Lb < 0.0 )
  g = @( x ) [ sum( x( x <= 0.0 ) ); sum( x( x >= 0.0 ) ) ] - [ Lb; Ub ];
else
  g = @( x ) sum( x ) - Ub;
endif

Everithing works, but ... one more thing - although initial guess is met equality constraints exactly, i.e. sum( R( R <= 0.0 ) ) = -1.0 and sum( R( R >= 0.0 ) ) = 1.0, the result of sqp solver is not - sum( x( x <= 0.0 ) ) is a bit less than -1.0 and sum( x( x >= 0.0 ) ) is a bit more than 1.0.
Do you have any ideas why it is so - maybe tune the tolerance of sqp?
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Newbie question - sqp solver equality constraints function

Juan Pablo Carbajal-2
On Thu, Jun 29, 2017 at 2:12 PM, alexegin <[hidden email]> wrote:
> Do you have any ideas why it is so - maybe tune the tolerance of sqp?

Yes, this is the tolerances, but I d not think you can change them
independently, there is a global tolerance parameter.
The output of sqp will give you information on how well the
constraints are respected, check the help.

Also, always try your ideas before asking. Do not abuse our time! :D

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Loading...