Mercurial > octave-nkf

changeset 6768:40e1255eda0e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-06-29 20:11:58 by jwe]
author jwe
date Fri, 29 Jun 2007 20:12:04 +0000
parents a6c8000f113e
children 77785733a18d
files scripts/ChangeLog scripts/optimization/sqp.m
diffstat 2 files changed, 117 insertions(+), 25 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ChangeLog	Thu Jun 28 19:42:42 2007 +0000
+++ b/scripts/ChangeLog	Fri Jun 29 20:12:04 2007 +0000
@@ -1,3 +1,9 @@
+2007-06-29  Marcus W. Reble  <reble@wisc.edu>
+
+	* optimization/sqp.m (sqp): New args, lb, ub, maxiter, and tolerance.
+	(fdjac): Set nx outside of if block.
+	(cf_ub_lb, cigrad_ub_lb): New subfunctons.
+
 2007-06-28  Michael Goffioul  <michael.goffioul@gmail.com>
 
 	* plot/subplot.m: Add 'ishandle' check when parsing the existing axes.
--- a/scripts/optimization/sqp.m	Thu Jun 28 19:42:42 2007 +0000
+++ b/scripts/optimization/sqp.m	Fri Jun 29 20:12:04 2007 +0000
@@ -182,14 +182,16 @@
 ## @seealso{qp}
 ## @end deftypefn
 
-function [x, obj, info, iter, nf, lambda] = sqp (x, objf, cef, cif)
+function [x, obj, info, iter, nf, lambda] = sqp (x, objf, cef, cif, lb, ub, maxiter, tolerance)
 
-  global nfun;
+  global __sqp_nfun__;
   global __sqp_obj_fun__;
   global __sqp_ce_fun__;
   global __sqp_ci_fun__;
+  global __sqp_cif__;
+  global __sqp_cifcn__;
 
-  if (nargin >= 2 && nargin <= 4)
+  if (nargin >= 2 && nargin <= 8 && nargin != 5)
 
     ## Choose an initial NxN symmetric positive definite Hessan
     ## approximation B.
@@ -245,29 +247,84 @@
 
     ci_fun = @empty_cf;
     ci_grd = @empty_jac;
+	
     if (nargin > 3)
-      ci_grd = @fd_ci_jac;
-      if (iscell (cif))
-	if (length (cif) > 0)
-	  __sqp_ci_fun__ = ci_fun = cif{1};
-	  if (length (cif) > 1)
-	    ci_grd = cif{2};
+      ## constraint function given by user with possibly gradient
+      __sqp_cif__ = cif;
+      ## constraint function given by user without gradient
+      __sqp_cifcn__ = @empty_cf;
+      if (iscell (__sqp_cif__))
+	if (length (__sqp_cif__) > 0)
+	  __sqp_cifcn__ = __sqp_cif__{1};
+	endif
+      elseif (! isempty (__sqp_cif__))
+	__sqp_cifcn__ = __sqp_cif__;
+      endif
+
+      if (nargin < 5)
+      	ci_grd = @fd_ci_jac;
+      	if (iscell (cif))
+	  if (length (cif) > 0)
+	    __sqp_ci_fun__ = ci_fun = cif{1};
+	    if (length (cif) > 1)
+	      ci_grd = cif{2};
+	    endif
+	  else
+	    error ("sqp: invalid equality constraint function");
 	  endif
+      	elseif (! isempty (cif))
+	  ci_fun = cif;
+      	endif
+      else
+	global __sqp_lb__;
+	if (isvector (lb))
+	  __sqp_lb__ = lb;
+	elseif (isempty (lb))
+	  __sqp_lb__ = -realmax;
 	else
-	  error ("sqp: invalid equality constraint function");
+	  error ("sqp: invalid lower bound");
 	endif
-      elseif (! isempty (cif))
-	ci_fun = cif;
+
+	global __sqp_ub__;
+	if (isvector (ub))
+	  __sqp_ub__ = ub;
+	elseif (isempty (lb))
+	  __sqp_ub__ = realmax;
+	else
+	  error ("sqp: invalid upper bound");
+	endif
+
+	if (lb > ub)
+	  error ("sqp: upper bound smaller than lower bound");
+	endif
+       	__sqp_ci_fun__ = ci_fun = @cf_ub_lb;
+       	ci_grd = @cigrad_ub_lb;
+      endif
+      __sqp_ci_fun__ = ci_fun;
+    endif
+
+    iter_max = 100;
+    if (nargin > 6 && ! isempty (maxiter))
+      if (isscalar (maxiter) && maxiter > 0 && round (maxiter) == maxiter))
+	iter_max = maxiter;
+      else
+	error ("sqp: invalid number of maximum iterations");
       endif
     endif
-    __sqp_ci_fun__ = ci_fun;
 
-    iter_max = 100;
+    tol = sqrt (eps);
+    if (nargin > 7 && ! isempty (tolerance))
+      if (isscalar (tolerance) && tolerance > 0)
+	tol = tolerance;
+      else
+	error ("sqp: invalid value for tolerance");
+      endif
+    endif
 
     iter = 0;
 
     obj = feval (obj_fun, x);
-    nfun = 1;
+    __sqp_nfun__ = 1;
 
     c = feval (obj_grd, x);
 
@@ -294,7 +351,7 @@
 
     ## report ();
 
-    ## report (iter, qp_iter, alpha, nfun, obj);
+    ## report (iter, qp_iter, alpha, __sqp_nfun__, obj);
 
     info = 0;
 
@@ -318,8 +375,6 @@
       t3 = all (lambda_i >= 0);
       t4 = norm (lambda .* con);
 
-      tol = sqrt (eps);
-
       if (t2 && t3 && max ([t0; t1; t4]) < tol)
 	break;
       endif
@@ -431,7 +486,7 @@
 
       A = A_new;
 
-      ## report (iter, qp_iter, alpha, nfun, obj);
+      ## report (iter, qp_iter, alpha, __sqp_nfun__, obj);
 
     endwhile
 
@@ -439,7 +494,7 @@
       info = 103;
     endif
 
-    nf = nfun;
+    nf = __sqp_nfun__;
 
   else
 
@@ -452,7 +507,7 @@
 
 function [merit, obj] = phi_L1 (obj, obj_fun, ce_fun, ci_fun, x, mu)
 
-  global nfun;
+  global __sqp_nfun__;
 
   ce = feval (ce_fun, x);
   ci = feval (ci_fun, x);
@@ -463,7 +518,7 @@
 
   if (isempty (obj))
     obj = feval (obj_fun, x);
-    nfun++;
+    __sqp_nfun__++;
   endif
 
   merit = obj;
@@ -564,7 +619,7 @@
 
 
 function jac = fdjac (f, x)
-
+  nx = length (x);
   if (! isempty (f))
     y0 = feval (f, x);
     nf = length (y0);
@@ -618,8 +673,39 @@
 
 function jac = fd_ci_jac (x)
 
-  global __sqp_ci_fun__;
+  global __sqp_cifcn__;
+  ## __sqp_cifcn__ = constraint function without gradients and lb or ub
+  jac = fdjac (__sqp_cifcn__, x);
+
+### endfunction
+
+function res = cf_ub_lb (x)
 
-  jac = fdjac (__sqp_ci_fun__, x);
+  ## combine constraint function with ub and lb
+  global __sqp_cifcn__ __sqp_lb__ __sqp_ub__
+
+  res = [x-__sqp_lb__; __sqp_ub__-x];
+
+  if (! isempty (__sqp_cifcn__))
+    res = [feval(__sqp_cifcn__,x); x-__sqp_lb__; __sqp_ub__-x];
+  endif
 
 ### endfunction
+
+function res = cigrad_ub_lb (x)
+
+  global __sqp_cif__
+
+  res = [eye(numel(x)); -eye(numel(x))];
+
+  cigradfcn = @fd_ci_jac;
+
+  if (iscell (__sqp_cif__) && length (__sqp_cif__) > 1)
+    cigradfcn = __sqp_cif__{2};
+  endif
+	
+  if (! isempty (cigradfcn))
+    res = [feval(cigradfcn,x); eye(numel(x)); -eye(numel(x))];
+  endif
+
+### endfunction
\ No newline at end of file