changelog shortlog tags changeset files revisions annotate raw

scripts/optimization/qp.m

changeset 10289: 4b124317dc38
parent:8f51a90eb8d1
author: John W. Eaton <jwe@octave.org>
date: Tue Feb 09 20:58:55 2010 -0500 (63 minutes ago)
permissions: -rw-r--r--
description: base_properties::set_children: account for hidden children 
     1## Copyright (C) 2000, 2001, 2004, 2005, 2006, 2007, 2008,

     2##               2009 Gabriele Pannocchia.

     3##

     4## This file is part of Octave.

     5##

     6## Octave is free software; you can redistribute it and/or modify it

     7## under the terms of the GNU General Public License as published by

     8## the Free Software Foundation; either version 3 of the License, or (at

     9## your option) any later version.

    10##

    11## Octave is distributed in the hope that it will be useful, but

    12## WITHOUT ANY WARRANTY; without even the implied warranty of

    13## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

    14## General Public License for more details.

    15##

    16## You should have received a copy of the GNU General Public License

    17## along with Octave; see the file COPYING.  If not, see

    18## <http://www.gnu.org/licenses/>.

    19

    20## -*- texinfo -*-

    21## @deftypefn {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@var{x0}, @var{H})

    22## @deftypefnx {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@var{x0}, @var{H}, @var{q})

    23## @deftypefnx {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@var{x0}, @var{H}, @var{q}, @var{A}, @var{b})

    24## @deftypefnx {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@var{x0}, @var{H}, @var{q}, @var{A}, @var{b}, @var{lb}, @var{ub})

    25## @deftypefnx {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@var{x0}, @var{H}, @var{q}, @var{A}, @var{b}, @var{lb}, @var{ub}, @var{A_lb}, @var{A_in}, @var{A_ub})

    26## @deftypefnx {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@dots{}, @var{options})

    27## Solve the quadratic program

    28## @tex

    29## $$

    30##  \min_x {1 \over 2} x^T H x + x^T q

    31## $$

    32## @end tex

    33## @ifnottex

    34##

    35## @example

    36## @group

    37##      min 0.5 x'*H*x + x'*q

    38##       x

    39## @end group

    40## @end example

    41##

    42## @end ifnottex

    43## subject to

    44## @tex

    45## $$

    46##  Ax = b \qquad lb \leq x \leq ub \qquad A_{lb} \leq A_{in} \leq A_{ub}

    47## $$

    48## @end tex

    49## @ifnottex

    50##

    51## @example

    52## @group

    53##      A*x = b

    54##      lb <= x <= ub

    55##      A_lb <= A_in*x <= A_ub

    56## @end group

    57## @end example

    58## @end ifnottex

    59##

    60## @noindent

    61## using a null-space active-set method.

    62##

    63## Any bound (@var{A}, @var{b}, @var{lb}, @var{ub}, @var{A_lb},

    64## @var{A_ub}) may be set to the empty matrix (@code{[]}) if not

    65## present.  If the initial guess is feasible the algorithm is faster.

    66##

    67## @table @var

    68## @item options

    69## An optional structure containing the following

    70## parameter(s) used to define the behavior of the solver.  Missing elements

    71## in the structure take on default values, so you only need to set the

    72## elements that you wish to change from the default.

    73##

    74## @table @code

    75## @item MaxIter (default: 200)

    76## Maximum number of iterations.

    77## @end table

    78## @end table

    79##

    80## The value @var{info} is a structure with the following fields:

    81## @table @code

    82## @item solveiter

    83## The number of iterations required to find the solution.

    84## @item info

    85## An integer indicating the status of the solution, as follows:

    86## @table @asis

    87## @item 0

    88## The problem is feasible and convex.  Global solution found.

    89## @item 1

    90## The problem is not convex.  Local solution found.

    91## @item 2

    92## The problem is not convex and unbounded.

    93## @item 3

    94## Maximum number of iterations reached.

    95## @item 6

    96## The problem is infeasible.

    97## @end table

    98## @end table

    99## @end deftypefn

   100

   101## PKG_ADD: __all_opts__ ("qp");

   102

   103function [x, obj, INFO, lambda] = qp (x0, H, varargin)

   104

   105  nargs = nargin;

   106

   107  if (nargin == 1 && ischar (x0) && strcmp (x0, 'defaults'))

   108    x = optimset ("MaxIter", 200);

   109    return;

   110  endif

   111

   112  if (nargs > 2 && isstruct (varargin{end}))

   113    options = varargin{end};

   114    nargs--;

   115  else

   116    options = struct ();

   117  endif

   118

   119  if (nargs >= 3)

   120    q = varargin{1};

   121  else

   122    q = [];

   123  endif

   124

   125  if (nargs >= 5)

   126    A = varargin{2};

   127    b = varargin{3};

   128  else

   129    A = [];

   130    b = [];

   131  endif

   132

   133  if (nargs >= 7)

   134    lb = varargin{4};

   135    ub = varargin{5};

   136  else

   137    lb = [];

   138    ub = [];

   139  endif

   140

   141  if (nargs == 10)

   142    A_lb = varargin{6};

   143    A_in = varargin{7};

   144    A_ub = varargin{8};

   145  else

   146    A_lb = [];

   147    A_in = [];

   148    A_ub = [];

   149  endif

   150

   151  if (nargs == 2 || nargs == 3 || nargs == 5 || nargs == 7 || nargs == 10)

   152

   153    maxit = optimget (options, "MaxIter", 200);

   154

   155    ## Checking the quadratic penalty

   156    if (! issquare (H))

   157      error ("qp: quadratic penalty matrix not square");

   158    elseif (! ishermitian (H))

   159      ## warning ("qp: quadratic penalty matrix not hermitian");

   160      H = (H + H')/2;

   161    endif

   162    n = rows (H);

   163

   164    ## Checking the initial guess (if empty it is resized to the

   165    ## right dimension and filled with 0)

   166    if (isempty (x0))

   167      x0 = zeros (n, 1);

   168    elseif (numel (x0) != n)

   169      error ("qp: the initial guess has incorrect length");

   170    endif

   171

   172    ## Linear penalty.

   173    if (isempty (q))

   174      q = zeros (n, 1);

   175    elseif (numel (q) != n)

   176      error ("qp: the linear term has incorrect length");

   177    endif

   178

   179    ## Equality constraint matrices

   180    if (isempty (A) || isempty (b))

   181      A = zeros (0, n);

   182      b = zeros (0, 1);

   183      n_eq = 0;

   184    else

   185      [n_eq, n1] = size (A);

   186      if (n1 != n)

   187	error ("qp: equality constraint matrix has incorrect column dimension");

   188      endif

   189      if (numel (b) != n_eq)

   190	error ("qp: equality constraint matrix and vector have inconsistent dimension");

   191      endif

   192    endif

   193

   194    ## Bound constraints

   195    Ain = zeros (0, n);

   196    bin = zeros (0, 1);

   197    n_in = 0;

   198    if (nargs > 5)

   199      if (! isempty (lb))

   200	if (numel (lb) != n)

   201	  error ("qp: lower bound has incorrect length");

   202	elseif (isempty (ub))

   203	  Ain = [Ain; eye(n)];

   204	  bin = [bin; lb];

   205	endif

   206      endif

   207

   208      if (! isempty (ub))

   209	if (numel (ub) != n)

   210	  error ("qp: upper bound has incorrect length");

   211	elseif (isempty (lb))

   212	  Ain = [Ain; -eye(n)];

   213	  bin = [bin; -ub];

   214	endif

   215      endif

   216      

   217      if (! isempty (lb) && ! isempty (ub))

   218	rtol = sqrt (eps);

   219	for i = 1:n

   220	  if (abs(lb (i) - ub(i)) < rtol*(1 + max (abs (lb(i) + ub(i)))))

   221            ## These are actually an equality constraint

   222	    tmprow = zeros(1,n);

   223	    tmprow(i) = 1;

   224            A = [A;tmprow];

   225            b = [b; 0.5*(lb(i) + ub(i))];

   226	    n_eq = n_eq + 1;

   227	  else

   228	    tmprow = zeros(1,n);

   229	    tmprow(i) = 1;

   230	    Ain = [Ain; tmprow; -tmprow];

   231	    bin = [bin; lb(i); -ub(i)];

   232	    n_in = n_in + 2;

   233	  endif

   234	endfor

   235      endif

   236    endif

   237

   238    ## Inequality constraints

   239    if (nargs > 7)

   240      [dimA_in, n1] = size (A_in);

   241      if (n1 != n)

   242	error ("qp: inequality constraint matrix has incorrect column dimension");

   243      else

   244	if (! isempty (A_lb))

   245	  if (numel (A_lb) != dimA_in)

   246	    error ("qp: inequality constraint matrix and lower bound vector inconsistent");

   247	  elseif (isempty (A_ub))

   248	    Ain = [Ain; A_in];

   249	    bin = [bin; A_lb];

   250	  endif

   251	endif

   252	if (! isempty (A_ub))

   253	  if (numel (A_ub) != dimA_in)

   254	    error ("qp: inequality constraint matrix and upper bound vector inconsistent");

   255	  elseif (isempty (A_lb))

   256	    Ain = [Ain; -A_in];

   257	    bin = [bin; -A_ub];

   258	  endif

   259	endif

   260	

   261	if (! isempty (A_lb) && ! isempty (A_ub))

   262	  rtol = sqrt (eps);

   263	  for i = 1:dimA_in

   264	    if (abs (A_lb(i) - A_ub(i)) < rtol*(1 + max (abs (A_lb(i) + A_ub(i)))))

   265              ## These are actually an equality constraint

   266	      tmprow = A_in(i,:);

   267              A = [A;tmprow];

   268              b = [b; 0.5*(A_lb(i) + A_ub(i))];

   269	      n_eq = n_eq + 1;

   270	    else

   271	      tmprow = A_in(i,:);

   272	      Ain = [Ain; tmprow; -tmprow];

   273	      bin = [bin; A_lb(i); -A_ub(i)];

   274	      n_in = n_in + 2;

   275	    endif

   276	  endfor

   277	endif

   278      endif

   279    endif

   280

   281    ## Now we should have the following QP:

   282    ##

   283    ##   min_x  0.5*x'*H*x + x'*q

   284    ##   s.t.   A*x = b

   285    ##          Ain*x >= bin

   286

   287    ## Discard inequality constraints that have -Inf bounds since those

   288    ## will never be active.

   289    idx = isinf (bin) & bin < 0;

   290

   291    bin(idx) = [];

   292    Ain(idx,:) = [];

   293

   294    n_in = numel (bin);

   295

   296    ## Check if the initial guess is feasible.

   297    if (isa (x0, "single") || isa (H, "single") || isa (q, "single") || isa (A, "single")

   298	|| isa (b, "single"))

   299      rtol = sqrt (eps ("single"));

   300    else

   301      rtol = sqrt (eps);

   302    endif

   303

   304    eq_infeasible = (n_eq > 0 && norm (A*x0-b) > rtol*(1+abs (b)));

   305    in_infeasible = (n_in > 0 && any (Ain*x0-bin < -rtol*(1+abs (bin))));

   306

   307    info = 0;

   308    if (eq_infeasible || in_infeasible)

   309      ## The initial guess is not feasible.

   310      ## First define xbar that is feasible with respect to the equality

   311      ## constraints.

   312      if (eq_infeasible)

   313	if (rank (A) < n_eq)

   314	  error ("qp: equality constraint matrix must be full row rank")

   315	endif

   316	xbar = pinv (A) * b;

   317      else

   318	xbar = x0;

   319      endif

   320

   321      ## Check if xbar is feasible with respect to the inequality

   322      ## constraints also.

   323      if (n_in > 0)

   324	res = Ain * xbar - bin;

   325	if (any (res < -rtol * (1 + abs (bin))))

   326	  ## xbar is not feasible with respect to the inequality

   327	  ## constraints.  Compute a step in the null space of the

   328	  ## equality constraints, by solving a QP.  If the slack is

   329	  ## small, we have a feasible initial guess.  Otherwise, the

   330	  ## problem is infeasible.

   331	  if (n_eq > 0)

   332	    Z = null (A);

   333	    if (isempty (Z))

   334	      ## The problem is infeasible because A is square and full

   335	      ## rank, but xbar is not feasible.

   336	      info = 6;

   337	    endif

   338	  endif

   339

   340	  if (info != 6)

   341            ## Solve an LP with additional slack variables to find

   342	    ## a feasible starting point.

   343	    gamma = eye (n_in);

   344	    if (n_eq > 0)

   345	      Atmp = [Ain*Z, gamma];

   346	      btmp = -res;

   347	    else

   348	      Atmp = [Ain, gamma];

   349	      btmp = bin;

   350	    endif

   351	    ctmp = [zeros(n-n_eq, 1); ones(n_in, 1)];

   352	    lb = [-Inf*ones(n-n_eq,1); zeros(n_in,1)];

   353	    ub = [];

   354	    ctype = repmat ("L", n_in, 1);

   355	    [P, dummy, status] = glpk (ctmp, Atmp, btmp, lb, ub, ctype);

   356	    if ((status == 180 || status == 181 || status == 151)

   357		&& all (abs (P(n-n_eq+1:end)) < rtol * (1 + norm (btmp))))

   358	      ## We found a feasible starting point

   359	      if (n_eq > 0)

   360		x0 = xbar + Z*P(1:n-n_eq);

   361	      else

   362		x0 = P(1:n);

   363              endif

   364	    else

   365	      ## The problem is infeasible

   366	      info = 6;

   367	    endif

   368	  endif

   369	else

   370	  ## xbar is feasible.  We use it a starting point.

   371	  x0 = xbar;

   372	endif

   373      else

   374	## xbar is feasible.  We use it a starting point.

   375	x0 = xbar;

   376      endif

   377    endif

   378

   379    if (info == 0)

   380      ## The initial (or computed) guess is feasible.

   381      ## We call the solver.

   382      [x, lambda, info, iter] = __qp__ (x0, H, q, A, b, Ain, bin, maxit);

   383    else

   384      iter = 0;

   385      x = x0;

   386      lambda = [];

   387    endif

   388    obj = 0.5 * x' * H * x + q' * x;

   389    INFO.solveiter = iter;

   390    INFO.info = info;

   391

   392  else

   393    print_usage ();

   394  endif

   395

   396endfunction