# HG changeset patch # User jwe # Date 1114094600 0 # Node ID b98bf1d70a0ae821f87fd52cfc11b5638468b751 # Parent 7434b60da83ee63dd0357769c6db74d229a16e67 [project @ 2005-04-21 14:43:14 by jwe] diff --git a/scripts/ChangeLog b/scripts/ChangeLog --- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,9 @@ +2005-04-21 John W. Eaton + + * optimization/glpk.m: Handle SENSE argument. + + * optimization/qp.m, optimization/sqp.m: New files. + 2005-04-08 John W. Eaton * Makefile.in (clean, distclean, maintainer-clean): diff --git a/scripts/optimization/glpk.m b/scripts/optimization/glpk.m --- a/scripts/optimization/glpk.m +++ b/scripts/optimization/glpk.m @@ -232,7 +232,7 @@ ## percentage of @code{tolbnd} or @code{toldj}). ## ## @item tolbnd (@code{LPX_K_TOLBND}, default: 10e-7) -## Relative tolerance used to check ifthe current basic solution is primal +## Relative tolerance used to check if the current basic solution is primal ## feasible. It is not recommended that you change this parameter unless you ## have a detailed understanding of its purpose. ## @@ -246,7 +246,7 @@ ## simplex table. It is not recommended that you change this parameter unless you ## have a detailed understanding of its purpose. ## -## @item objll (@code{LPX_K_OBJLL)}, default: -DBL_MAX +## @item objll (@code{LPX_K_OBJLL}, default: -DBL_MAX) ## Lower limit of the objective function. If on the phase II the objective ## function reaches this limit and continues decreasing, the solver stops ## the search. This parameter is used in the dual simplex method only. @@ -268,7 +268,7 @@ ## output. Non-positive value means no delay. ## ## @item tolint (@code{LPX_K_TOLINT}, default: 10e-5) -## Relative tolerance used to check ifthe current basic solution is integer +## Relative tolerance used to check if the current basic solution is integer ## feasible. It is not recommended that you change this parameter unless ## you have a detailed understanding of its purpose. ## @@ -477,7 +477,7 @@ ## 7) Vector with the type of variables if (nargin > 6) - if isempty (vartype) + if (isempty (vartype)) vartype = repmat ("C", nx, 1); elseif (! ischar (vartype) || all (size (vartype) > 1) || length (vartype) != nx) @@ -492,7 +492,23 @@ vartype = repmat ("C", nx, 1); endif - ## 8) Parameters vector + ## 8) Sense of optimization + + if (nargin > 7) + if (isempty (sense)) + sense = 1; + elseif (ischar (sense) || all (size (sense) > 1) || ! isreal (sense)) + error ("SENSE must be an integer value"); + elseif (sense >= 0) + sense = 1; + else + sense = -1; + endif + else + sense = 1; + endif + + ## 9) Parameters vector if (nargin > 8) if (! isstruct (param)) diff --git a/scripts/optimization/qp.m b/scripts/optimization/qp.m new file mode 100644 --- /dev/null +++ b/scripts/optimization/qp.m @@ -0,0 +1,280 @@ +## Copyright (C) 2000, 2001, 2004, 2005 Gabriele Pannocchia. +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, 59 Temple Place - Suite 330, Boston, MA +## 02111-1307, USA. + +## -*- texinfo -*- +## @deftypefn {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{Ain}, @var{A_ub}) +## Solve the quadratic program +## @ifinfo +## +## @example +## min 0.5 x'*H*x + x'*q +## x +## @end example +## +## @end ifinfo +## @iftex +## @tex +## @end tex +## @end iftex +## subject to +## @ifinfo +## +## @example +## A*x = b +## lb <= x <= ub +## A_lb <= Ain*x <= A_ub +## @end example +## @end ifinfo +## @iftex +## @tex +## @end tex +## @end iftex +## +## @noindent +## using a null-space active-set method. +## +## Any bound (@var{A}, @var{b}, @var{lb}, @var{ub}, @var{A_lb}, +## @var{A_ub}) may be set to the empty matrix (@code{[]}) if not +## present. If the initial guess is feasible the algorithm is faster. +## +## The value @var{info} is a structure with the following fields: +## @table @code +## @item solveiter +## The number of iterations required to find the solution. +## @item info +## An integer indicating the status of the solution, as follows: +## @table @asis +## @item 0 +## The problem is feasible and convex. Global solution found. +## @item 1 +## The problem is not convex. Local solution found. +## @item 2 +## The problem is not convex and unbounded. +## @item 3 +## Maximum number of iterations reached. +## @item 6 +## The problem is infeasible. +## @end table +## @end table +## @end deftypefn + +function [x, obj, INFO, lambda] = qp (x0, H, q, A, b, lb, ub, A_lb, A_in, A_ub) + + if (nargin == 5 || nargin == 7 || nargin == 10) + + ## Checking the quadratic penalty + n = issquare (H); + if (n == 0) + error ("qp: quadratic penalty matrix not square"); + endif + + n1 = issymmetric (H); + if (n1 == 0) + ## warning ("qp: quadratic penalty matrix not symmetric"); + H = (H + H')/2; + endif + + ## Checking the initial guess (if empty it is resized to the + ## right dimension and filled with 0) + if (isempty (x0)) + x0 = zeros (n, 1); + elseif (length (x0) != n) + error ("qp: the initial guess has incorrect length"); + endif + + ## Linear penalty. + if (length (q) != n) + error ("qp: the linear term has incorrect length"); + endif + + ## Equality constraint matrices + if (isempty (A) || isempty(b)) + n_eq = 0; + A = zeros (n_eq, n); + b = zeros (n_eq, 1); + else + [n_eq, n1] = size (A); + if (n1 != n) + error ("qp: equality constraint matrix has incorrect column dimension"); + endif + if (length (b) != n_eq) + error ("qp: equality constraint matrix and vector have inconsistent dimension"); + endif + endif + + ## Bound constraints + Ain = zeros (0, n); + bin = zeros (0, 1); + n_in = 0; + if (nargin > 5) + if (! isempty (lb)) + if (length(lb) != n) + error ("qp: lower bound has incorrect length"); + else + Ain = [Ain; eye(n)]; + bin = [bin; lb]; + endif + endif + + if (! isempty (ub)) + if (length (ub) != n) + error ("qp: upper bound has incorrect length"); + else + Ain = [Ain; -eye(n)]; + bin = [bin; -ub]; + endif + endif + endif + + ## Inequality constraints + if (nargin > 7) + [dimA_in, n1] = size (A_in); + if (n1 != n) + error ("qp: inequality constraint matrix has incorrect column dimension"); + else + if (! isempty (A_lb)) + if (length (A_lb) != dimA_in) + error ("qp: inequality constraint matrix and lower bound vector inconsistent"); + else + Ain = [Ain; A_in]; + bin = [bin; A_lb]; + endif + endif + if (! isempty (A_ub)) + if (length (A_ub) != dimA_in) + error ("qp: inequality constraint matrix and upper bound vector inconsistent"); + else + Ain = [Ain; -A_in]; + bin = [bin; -A_ub]; + endif + endif + endif + endif + n_in = length (bin); + + ## Now we should have the following QP: + ## + ## min_x 0.5*x'*H*x + x'*q + ## s.t. A*x = b + ## Ain*x >= bin + + ## Discard inequality constraints that have -Inf bounds since those + ## will never be active. + idx = isinf (bin) & bin < 0; + bin(idx) = []; + Ain(idx,:) = []; + + ## Check if the initial guess is feasible. + rtol = sqrt (eps); + + eq_infeasible = (n_eq > 0 && norm (A*x0-b) > rtol*(1+norm (b))); + in_infeasible = (n_in > 0 && any (Ain*x0-bin < -rtol*(1+norm (bin)))); + + info = 0; + if (eq_infeasible || in_infeasible) + ## The initial guess is not feasible. + ## First define xbar that is feasible with respect to the equality + ## constraints. + if (eq_infeasible) + if (rank (A) < n_eq) + error ("qp: equality constraint matrix must be full row rank") + endif + xbar = pinv (A) * b; + else + xbar = x0; + endif + + ## Check if xbar is feasible with respect to the inequality + ## constraints also. + if (n_in > 0) + res = Ain * xbar - bin; + if (any (res < -rtol * (1 + norm (bin)))) + ## xbar is not feasible with respect to the inequality + ## constraints. Compute a step in the null space of the + ## equality constraints, by solving a QP. If the slack is + ## small, we have a feasible initial guess. Otherwise, the + ## problem is infeasible. + if (n_eq > 0) + Z = null (A); + if (isempty (Z)) + ## The problem is infeasible because A is square and full + ## rank, but xbar is not feasible. + info = 6; + endif + endif + + if (info != 6) + ## Solve an LP with additional slack variables to find + ## a feasible starting point. + gamma = eye (n_in); + if (n_eq > 0) + Atmp = [Ain*Z, gamma]; + btmp = -res; + else + Atmp = [Ain, gamma]; + btmp = bin; + endif + ctmp = [zeros(n-n_eq, 1); ones(n_in, 1)]; + lb = [-Inf*ones(n-n_eq,1); zeros(n_in,1)]; + ub = []; + ctype = repmat ("L", n_in, 1); + [P, dummy, status] = glpk (ctmp, Atmp, btmp, lb, ub, ctype); + + if ((status == 180 || status == 181 || status == 151) + && all (abs (P(n-n_eq+1:end)) < rtol * (1 + norm (btmp)))) + ## We found a feasible starting point + if (n_eq > 0) + x0 = xbar + Z*P(1:n-n_eq) + else + x0 = P(1:n); + endif + else + ## The problem is infeasible + info = 6; + endif + endif + else + ## xbar is feasible. We use it a starting point. + x0 = xbar; + endif + else + ## xbar is feasible. We use it a starting point. + x0 = xbar; + endif + endif + + if (info == 0) + ## The initial (or computed) guess is feasible. + ## We call the solver. + maxit = 100; + [x, lambda, info, iter] = __qp__ (x0, H, q, A, b, Ain, bin, maxit); + else + iter = 0; + x = x0; + lambda = []; + endif + obj = 0.5 * x' * H * x + q' * x; + INFO.solveiter = iter; + INFO.info = info; + + else + usage ("[x, obj, info, lambda] = qp (x0, H, q, A, b, lb, ub, A_lb, Ain, A_ub)"); + endif + +endfunction diff --git a/scripts/optimization/sqp.m b/scripts/optimization/sqp.m new file mode 100644 --- /dev/null +++ b/scripts/optimization/sqp.m @@ -0,0 +1,603 @@ +## Copyright (C) 2005 John W. Eaton +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, 59 Temple Place - Suite 330, Boston, MA +## 02111-1307, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{x}, @var{obj}, @var{info}, @var{iter}, @var{nf}, @var{lambda}] =} sqp (@var{x}, @var{phi}, @var{g}, @var{h}) +## Solve the nonlinear program +## @ifinfo +## +## @example +## min phi (x) +## x +## @end example +## +## @end ifinfo +## @iftex +## @tex +## @end tex +## @end iftex +## subject to +## @ifinfo +## +## @example +## g(x) = 0 +## h(x) >= 0 +## @end example +## @end ifinfo +## @iftex +## @tex +## @end tex +## @end iftex +## +## @noindent +## using a successive quadratic programming method. +## +## The first argument is the initial guess for the vector @var{x}. +## +## The second argument is a function handle pointing to the ojective +## function. The objective function must be of the form +## +## @example +## y = phi (x) +## @end example +## +## @noindent +## in which @var{x} is a vector and @var{y} is a scalar. +## +## The second argument may also be a 2- or 3-element cell array of +## function handles. The first element should point to the objective +## function, the second should point to a function that computes the +## gradient of the objective function, and the third should point to a +## function to compute the hessian of the objective function. If the +## gradient function is not supplied, the gradient is computed by finite +## differences. If the hessian function is not supplied, a BFGS update +## formula is used to approximate the hessian. +## +## If supplied, the gradient function must be of the form +## +## @example +## g = gradient (x) +## @end example +## +## @noindent +## in which @var{x} is a vector and @var{g} is a vector. +## +## If supplied, the hessian function must be of the form +## +## @example +## h = hessian (x) +## @end example +## +## @noindent +## in which @var{x} is a vector and @var{h} is a matrix. +## +## The third and fourth arguments are function handles pointing to +## functions that compute the equality constraints and the inequality +## constraints, respectively. +## +## If your problem does not have equality (or inequality) constraints, +## you may pass an empty matrix for @var{cef} (or @var{cif}). +## +## If supplied, the equality and inequality constraint functions must be +## of the form +## +## @example +## r = f (x) +## @end example +## +## @noindent +## in which @var{x} is a vector and @var{r} is a vector. +## +## The third and fourth arguments may also be 2-element cell arrays of +## function handles. The first element should point to the constraint +## function and the second should point to a function that computes the +## gradient of the constraint function: +## +## @example +## [ d f(x) d f(x) d f(x) ] +## transpose ( [ ------ ----- ... ------ ] ) +## [ dx_1 dx_2 dx_N ] +## @end example +## +## Here is an example of calling @code{sqp}: +## +## @example +## function r = g (x) +## r = [ sumsq(x)-10; x(2)*x(3)-5*x(4)*x(5); x(1)^3+x(2)^3+1]; +## endfunction +## +## function obj = phi (x) +## obj = exp(prod(x)) - 0.5*(x(1)^3+x(2)^3+1)^2; +## endfunction +## +## x0 = [-1.8; 1.7; 1.9; -0.8; -0.8]; +## +## [x, obj, info, iter, nf, lambda] = sqp (x0, @@phi, @@g, []) +## +## x = +## +## -1.71714 +## 1.59571 +## 1.82725 +## -0.76364 +## -0.76364 +## +## obj = 0.053950 +## info = 101 +## iter = 8 +## nf = 10 +## lambda = +## +## -0.0401627 +## 0.0379578 +## -0.0052227 +## @end example +## +## The value returned in @var{info} may be one of the following: +## @table @asis +## @item 101 +## The algorithm terminated because the norm of the last step was less +## than @code{tol * norm (x))} (the value of tol is currently fixed at +## @code{sqrt (eps)}---edit @file{sqp.m} to modify this value. +## @item 102 +## The BFGS update failed. +## @item 103 +## The maximum number of iterations was reached (the maximum number of +## allowed iterations is currently fixed at 100---edit @file{sqp.m} to +## increase this value). +## @end table +## @end deftypefn +## @seealso{qp} + +function [x, obj, info, iter, nf, lambda] = sqp (x, objf, cef, cif) + + global nfun; + global __sqp_obj_fun__; + global __sqp_ce_fun__; + global __sqp_ci_fun__; + + if (nargin >= 2 && nargin <= 4) + + ## Choose an initial NxN symmetric positive definite Hessan + ## approximation B. + + n = length (x); + + B = eye (n, n); + + ## Evaluate objective function, constraints, and gradients at initial + ## value of x. + ## + ## obj_fun + ## obj_grad + ## ce_fun -- equality constraint functions + ## ci_fun -- inequality constraint functions + ## A == [grad_{x_1} cx_fun, grad_{x_2} cx_fun, ..., grad_{x_n} cx_fun]^T + + obj_grd = @fd_obj_grd; + have_hess = 0; + if (iscell (objf)) + if (length (objf) > 0) + __sqp_obj_fun__ = obj_fun = objf{1}; + if (length (objf) > 1) + obj_grd = objf{2}; + if (length (objf) > 2) + obj_hess = objf{3}; + have_hess = 1; + B = feval (obj_hess, x); + endif + endif + else + error ("sqp: invalid objective function"); + endif + else + __sqp_obj_fun__ = obj_fun = objf; + endif + + ce_fun = @empty_cf; + ce_grd = @empty_jac; + if (nargin > 2) + ce_grd = @fd_ce_jac; + if (iscell (cef)) + if (length (cef) > 0) + __sqp_ce_fun__ = ce_fun = cef{1}; + if (length (cef) > 1) + ce_grd = cef{2}; + endif + else + error ("sqp: invalid equality constraint function"); + endif + elseif (! isempty (cef)) + ce_fun = cef; + endif + endif + __sqp_ce_fun__ = ce_fun; + + 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}; + endif + else + error ("sqp: invalid equality constraint function"); + endif + elseif (! isempty (cif)) + ci_fun = cif; + endif + endif + __sqp_ci_fun__ = ci_fun; + + iter_max = 100; + + iter = 0; + + obj = feval (obj_fun, x); + nfun = 1; + + c = feval (obj_grd, x); + + ce = feval (ce_fun, x); + F = feval (ce_grd, x); + + ci = feval (ci_fun, x); + C = feval (ci_grd, x); + + A = [F; C]; + + ## Choose an initial lambda (x is provided by the caller). + + lambda = 100 * ones (rows (A), 1); + + qp_iter = 1; + alpha = 1; + + ## report (); + + ## report (iter, qp_iter, alpha, nfun, obj); + + while (++iter < iter_max) + + ## Check convergence. This is just a simple check on the first + ## order necessary conditions. + + ## IDX is the indices of the active inequality constraints. + + nr_f = rows (F); + + lambda_e = lambda((1:nr_f)'); + lambda_i = lambda((nr_f+1:end)'); + + con = [ce; ci]; + + t0 = norm (c - A' * lambda); + t1 = norm (ce); + t2 = all (ci >= 0); + t3 = all (lambda_i >= 0); + t4 = norm (lambda .* con); + + tol = sqrt (eps); + + if (t2 && t3 && max ([t0; t1; t4]) < tol) + break; + endif + + ## Compute search direction p by solving QP. + + g = -ce; + d = -ci; + + ## Discard inequality constraints that have -Inf bounds since those + ## will never be active. + idx = isinf (d) & d < 0; + d(idx) = []; + C(idx,:) = []; + + [p, obj_qp, INFO, lambda] = qp (x, B, c, F, g, [], [], d, C, + Inf * ones (size (d))); + + info = INFO.info; + + ## Check QP solution and attempt to recover if it has failed. + + ## Choose mu such that p is a descent direction for the chosen + ## merit function phi. + + [x_new, alpha, obj_new] = linesearch_L1 (x, p, obj_fun, obj_grd, + ce_fun, ci_fun, lambda, obj); + + ## Evaluate objective function, constraints, and gradients at + ## x_new. + + c_new = feval (obj_grd, x_new); + + ce_new = feval (ce_fun, x_new); + F_new = feval (ce_grd, x_new); + + ci_new = feval (ci_fun, x_new); + C_new = feval (ci_grd, x_new); + + A_new = [F_new; C_new]; + + ## Set + ## + ## s = alpha * p + ## y = grad_x L (x_new, lambda) - grad_x L (x, lambda}) + + y = c_new - c; + + if (! isempty (A)) + t = ((A_new - A)'*lambda); + y -= t; + endif + + delx = x_new - x; + + if (norm (delx) < tol * norm (x)) + info = 101; + break; + endif + + if (have_hess) + + B = feval (obj_hess, x); + + else + + ## Update B using a quasi-Newton formula. + + delxt = delx'; + + ## Damped BFGS. Or maybe we would actually want to use the Hessian + ## of the Lagrangian, computed directly. + + d1 = delxt*B*delx; + + t1 = 0.2 * d1; + t2 = delxt*y; + + if (t2 < t1) + theta = 0.8*d1/(d1 - t2); + else + theta = 1; + endif + + r = theta*y + (1-theta)*B*delx; + + d2 = delxt*r; + + if (d1 == 0 || d2 == 0) + info = 102; + break; + endif + + B = B - B*delx*delxt*B/d1 + r*r'/d2; + + endif + + x = x_new; + + obj = obj_new; + + c = c_new; + + ce = ce_new; + F = F_new; + + ci = ci_new; + C = C_new; + + A = A_new; + + ## report (iter, qp_iter, alpha, nfun, obj); + + endwhile + + if (iter >= iter_max) + info = 103; + endif + + nf = nfun; + + else + + usage ("[x, obj, info, iter, nf, lambda] = sqp (x, objf, cef, cif)"); + + endif + +### endfunction + + +function [merit, obj] = phi_L1 (obj, obj_fun, ce_fun, ci_fun, x, mu) + + global nfun; + + ce = feval (ce_fun, x); + ci = feval (ci_fun, x); + + idx = ci < 0; + + con = [ce; ci(idx)]; + + if (isempty (obj)) + obj = feval (obj_fun, x); + nfun++; + endif + + merit = obj; + t = norm (con, 1) / mu; + + if (! isempty (t)) + merit += t; + endif + +### endfunction + + +function [x_new, alpha, obj] = linesearch_L1 (x, p, obj_fun, obj_grd, + ce_fun, ci_fun, lambda, obj) + + ## Choose parameters + ## + ## eta in the range (0, 0.5) + ## tau in the range (0, 1) + + eta = 0.25; + tau = 0.5; + + delta_bar = sqrt (eps); + + if (isempty (lambda)) + mu = 1 / delta_bar; + else + mu = 1 / (norm (lambda, Inf) + delta_bar); + endif + + alpha = 1; + + c = feval (obj_grd, x); + ce = feval (ce_fun, x); + + [phi_x_mu, obj] = phi_L1 (obj, obj_fun, ce_fun, ci_fun, x, mu); + + D_phi_x_mu = c' * p; + d = feval (ci_fun, x); + ## only those elements of d corresponding + ## to violated constraints should be included. + idx = d < 0; + t = - norm ([ce; d(idx)], 1) / mu; + if (! isempty (t)) + D_phi_x_mu += t; + endif + + while (1) + [p1, obj] = phi_L1 ([], obj_fun, ce_fun, ci_fun, x+alpha*p, mu); + p2 = phi_x_mu+eta*alpha*D_phi_x_mu; + if (p1 > p2) + ## Reset alpha = tau_alpha * alpha for some tau_alpha in the + ## range (0, tau). + tau_alpha = 0.9 * tau; ## ?? + alpha = tau_alpha * alpha; + else + break; + endif + endwhile + + ## Set x_new = x + alpha * p; + + x_new = x + alpha * p; + +### endfunction + + +function report (iter, qp_iter, alpha, nfun, obj) + + if (nargin == 0) + printf (" Itn ItQP Step Nfun Objective\n"); + else + printf ("%5d %4d %8.1g %5d %13.6e\n", iter, qp_iter, alpha, nfun, obj); + endif + +### endfunction + + +function grd = fdgrd (f, x) + + if (! isempty (f)) + y0 = feval (f, x); + nx = length (x); + grd = zeros (nx, 1); + deltax = sqrt (eps); + for i = 1:nx + t = x(i); + x(i) += deltax; + grd(i) = (feval (f, x) - y0) / deltax; + x(i) = t; + endfor + else + grd = zeros (0, 1); + endif + +### endfunction + + +function jac = fdjac (f, x) + + if (! isempty (f)) + y0 = feval (f, x); + nf = length (y0); + nx = length (x); + jac = zeros (nf, nx); + deltax = sqrt (eps); + for i = 1:nx + t = x(i); + x(i) += deltax; + jac(:,i) = (feval (f, x) - y0) / deltax; + x(i) = t; + endfor + else + jac = zeros (0, nx); + endif + +### endfunction + + +function grd = fd_obj_grd (x) + + global __sqp_obj_fun__; + + grd = fdgrd (__sqp_obj_fun__, x); + +### endfunction + + +function res = empty_cf (x) + + res = zeros (0, 1); + +### endfunction + + +function res = empty_jac (x) + + res = zeros (0, length (x)); + +### endfunction + + +function jac = fd_ce_jac (x) + + global __sqp_ce_fun__; + + jac = fdjac (__sqp_ce_fun__, x); + +### endfunction + + +function jac = fd_ci_jac (x) + + global __sqp_ci_fun__; + + jac = fdjac (__sqp_ci_fun__, x); + +### endfunction