fixing error status for PATH hook
Michel Juillard
parent
11161063d4
commit
450d7b099a
|
|
@ -206,7 +206,8 @@ elseif options.solve_algo == 11
|
|||
end
|
||||
olcppath = options.lcppath;
|
||||
[junk,M] = func(x,varargin{:});
|
||||
[x,mu] = pathlcp(fjac,olcppath.q,olcppath.lb,olcppath.ub,x,olcppath.A,olcppath.b,olcppath.t,olcppath.mu0);
|
||||
[x,mu,status] = pathlcp(fjac,olcppath.q,olcppath.lb,olcppath.ub,x,olcppath.A,olcppath.b,olcppath.t,olcppath.mu0);
|
||||
info = ~status;
|
||||
elseif options.solve_algo == 12
|
||||
% PATH mixed complementary problem
|
||||
% PATH linear mixed complementary problem
|
||||
|
|
@ -219,7 +220,8 @@ elseif options.solve_algo == 12
|
|||
global mcp_data
|
||||
mcp_data.func = func;
|
||||
mcp_data.args = varargin;
|
||||
[x,mu] = pathmcp(x,omcppath.lb,omcppath.ub,'mcp_func',omcppath.A,omcppath.b,omcppath.t,omcppath.mu0);
|
||||
[x,fval,jac,mu,status] = pathmcp(x,omcppath.lb,omcppath.ub,'mcp_func',omcppath.A,omcppath.b,omcppath.t,omcppath.mu0);
|
||||
info = ~status;
|
||||
else
|
||||
error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10:12]')
|
||||
end
|
||||
|
|
|
|||
|
|
@ -0,0 +1,172 @@
|
|||
function [z,mu,status] = pathlcp(M,q,l,u,z,A,b,t,mu)
|
||||
% pathlcp(M,q,l,u,z,A,b,t,mu)
|
||||
%
|
||||
% Solve the standard linear complementarity problem using PATH:
|
||||
% z >= 0, Mz + q >= 0, z'*(Mz + q) = 0
|
||||
%
|
||||
% Required input:
|
||||
% M(n,n) - matrix
|
||||
% q(n) - vector
|
||||
%
|
||||
% Output:
|
||||
% z(n) - solution
|
||||
% mu(m) - multipliers (if polyhedral constraints are present)
|
||||
%
|
||||
% Optional input:
|
||||
% l(n) - lower bounds default: zero
|
||||
% u(n) - upper bounds default: infinity
|
||||
% z(n) - starting point default: zero
|
||||
% A(m,n) - polyhedral constraint matrix default: empty
|
||||
% b(m) - polyhedral right-hand side default: empty
|
||||
% t(m) - type of polyhedral constraint default: 1
|
||||
% < 0: less than or equal
|
||||
% 0: equation
|
||||
% > 0: greater than or equal
|
||||
% mu(m) - starting value for multipliers default: zero
|
||||
%
|
||||
% The optional lower and upper bounds are used to define a linear mixed
|
||||
% complementarity problem (box constrained variational inequality).
|
||||
% l <= z <= u
|
||||
% where l_i < z_i < u_i => (Mz + q)_i = 0
|
||||
% l_i = z => (Mz + q)_i >= 0
|
||||
% u_i = z => (Mz + q)_i <= 0
|
||||
%
|
||||
% The optional constraints are used to define a polyhedrally constrained
|
||||
% variational inequality. These are transformed internally to a standard
|
||||
% mixed complementarity problem. The polyhedral constraints are of the
|
||||
% form
|
||||
% Ax ? b
|
||||
% where ? can be <=, =, or >= depending on the type specified for each
|
||||
% constraint.
|
||||
|
||||
Big = 1e20;
|
||||
|
||||
if (nargin < 2)
|
||||
error('two input arguments required for lcp(M, q)');
|
||||
end
|
||||
|
||||
if (~issparse(M))
|
||||
M = sparse(M); % Make sure M is sparse
|
||||
end
|
||||
q = full(q(:)); % Make sure q is a column vector
|
||||
|
||||
[mm,mn] = size(M); % Get the size of the inputs
|
||||
n = length(q);
|
||||
|
||||
if (mm ~= mn | mm ~= n)
|
||||
error('dimensions of M and q must match');
|
||||
end
|
||||
|
||||
if (n == 0)
|
||||
error('empty model');
|
||||
end
|
||||
|
||||
if (nargin < 3 | isempty(l))
|
||||
l = zeros(n,1);
|
||||
end
|
||||
|
||||
if (nargin < 4 | isempty(u))
|
||||
u = Big*ones(n,1);
|
||||
end
|
||||
|
||||
if (nargin < 5 | isempty(z))
|
||||
z = zeros(n,1);
|
||||
end
|
||||
|
||||
z = full(z(:)); l = full(l(:)); u = full(u(:));
|
||||
if (length(z) ~= n | length(l) ~= n | length(u) ~= n)
|
||||
error('Input arguments are of incompatible sizes');
|
||||
end
|
||||
|
||||
l = max(l,-Big*ones(n,1));
|
||||
u = min(u,Big*ones(n,1));
|
||||
z = min(max(z,l),u);
|
||||
|
||||
m = 0;
|
||||
if (nargin > 5)
|
||||
if (nargin < 7)
|
||||
error('Polyhedral constraints require A and b');
|
||||
end
|
||||
|
||||
if (~issparse(A))
|
||||
A = sparse(A);
|
||||
end
|
||||
b = full(b(:));
|
||||
|
||||
m = length(b);
|
||||
|
||||
if (m > 0)
|
||||
|
||||
[am, an] = size(A);
|
||||
|
||||
if (am ~= m | an ~= n)
|
||||
error('Polyhedral constraints of incompatible sizes');
|
||||
end
|
||||
|
||||
if (nargin < 8 | isempty(t))
|
||||
t = ones(m,1);
|
||||
end
|
||||
|
||||
if (nargin < 9 | isempty(mu))
|
||||
mu = zeros(m,1);
|
||||
end
|
||||
|
||||
t = full(t(:)); mu = full(mu(:));
|
||||
if (length(t) ~= m | length(mu) ~= m)
|
||||
error('Polyhedral input arguments are of incompatible sizes');
|
||||
end
|
||||
|
||||
l_p = -Big*ones(m,1);
|
||||
u_p = Big*ones(m,1);
|
||||
|
||||
idx = find(t > 0);
|
||||
if (length(idx) > 0)
|
||||
l_p(idx) = zeros(length(idx),1);
|
||||
end
|
||||
|
||||
idx = find(t < 0);
|
||||
if (length(idx) > 0)
|
||||
u_p(idx) = zeros(length(idx),1);
|
||||
end
|
||||
|
||||
mu = min(max(mu,l_p),u_p);
|
||||
|
||||
M = [M -A'; A sparse(m,m)];
|
||||
q = [q; -b];
|
||||
|
||||
z = [z; mu];
|
||||
l = [l; l_p];
|
||||
u = [u; u_p];
|
||||
else
|
||||
if (nargin >= 9 & ~isempty(mu))
|
||||
error('No polyhedral constraints -- multipliers set.');
|
||||
end
|
||||
|
||||
if (nargin >= 8 & ~isempty(t))
|
||||
error('No polyhedral constraints -- equation types set.');
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
idx = find(l > u);
|
||||
if length(idx) > 0
|
||||
error('Bounds infeasible.');
|
||||
end
|
||||
|
||||
nnzJ = nnz(M);
|
||||
|
||||
[status, ttime] = lcppath(n+m, nnzJ, z, l, u, M, q);
|
||||
|
||||
%if (status ~= 1)
|
||||
% status,
|
||||
% error('Path fails to solve problem');
|
||||
%end
|
||||
|
||||
mu = [];
|
||||
if (m > 0)
|
||||
mu = z(n+1:n+m);
|
||||
z = z(1:n);
|
||||
end
|
||||
|
||||
return;
|
||||
|
||||
|
|
@ -0,0 +1,197 @@
|
|||
function [z,f,J,mu,status] = pathmcp(z,l,u,cpfj,A,b,t,mu)
|
||||
% pathmcp(z,l,u,cpfj,A,b,t,mu)
|
||||
%
|
||||
% Solve a polyhedrally constrained variational inequality using PATH
|
||||
%
|
||||
% Calling syntax: [z,f,J] = pathmcp(z,l,u,cpfunjac,A,b,t,mu)
|
||||
%
|
||||
% Input:
|
||||
% z - starting point
|
||||
% l - lower bounds on z
|
||||
% u - upper bounds on z
|
||||
%
|
||||
% cpfunjac - the name of the m-file for evaluating the function F and its
|
||||
% Jacobian J (without .m-extension).
|
||||
%
|
||||
% The following m-file must be supplied (where default name is
|
||||
% 'mcp_funjac.m' unless stated otherwise in the variable cpfunjac).
|
||||
%
|
||||
% 'mcp_funjac.m' contains function [f,J,domerr]=cpfunjac(z,jacflag)
|
||||
% that computes the function F and if jacflag=1 the sparse
|
||||
% Jacobian J at the point z. domerr returns the number of domain
|
||||
% violations.
|
||||
%
|
||||
% A - constraint matrix
|
||||
% b - right hand side of the constraints
|
||||
% t - types of the constraints
|
||||
% <0 : less than or equal
|
||||
% =0 : equal to
|
||||
% >0 : greater than or equal
|
||||
%
|
||||
% We have Ax ? b, ? is the type of constraint
|
||||
%
|
||||
% Output:
|
||||
% z - solution
|
||||
% mu - multipliers on the constraints
|
||||
% f - function evaluation at the solution
|
||||
% J - jacobian evaluation at the solution
|
||||
|
||||
|
||||
Big = 1e20;
|
||||
|
||||
if (nargin < 1)
|
||||
error('one input arguments required for mcp(z)');
|
||||
end
|
||||
|
||||
z = full(z(:));
|
||||
n = length(z);
|
||||
|
||||
if (n == 0)
|
||||
error('empty model');
|
||||
end
|
||||
|
||||
if (nargin < 2 | isempty(l))
|
||||
l = zeros(n,1);
|
||||
end
|
||||
|
||||
if (nargin < 3 | isempty(u))
|
||||
u = Big*ones(n,1);
|
||||
end
|
||||
|
||||
l = full(l(:)); u = full(u(:));
|
||||
if (length(l) ~= n | length(u) ~= n)
|
||||
error('Input arguments are of incompatible sizes');
|
||||
end
|
||||
|
||||
l = max(l,-Big*ones(n,1));
|
||||
u = min(u,Big*ones(n,1));
|
||||
z = min(max(z,l),u);
|
||||
|
||||
if (nargin < 4 | isempty(cpfj))
|
||||
cpfj = 'mcp_funjac';
|
||||
end
|
||||
|
||||
m = 0;
|
||||
mu = [];
|
||||
l_p = [];
|
||||
u_p = [];
|
||||
|
||||
if (nargin > 4)
|
||||
if (nargin < 6)
|
||||
error('Polyhedral constraints require A and b');
|
||||
end
|
||||
|
||||
if (~issparse(A))
|
||||
A = sparse(A);
|
||||
end
|
||||
b = full(b(:));
|
||||
|
||||
m = length(b);
|
||||
|
||||
if (m > 0)
|
||||
|
||||
[am, an] = size(A);
|
||||
|
||||
if (am ~= m | an ~= n)
|
||||
error('Polyhedral constraints of incompatible sizes');
|
||||
end
|
||||
|
||||
if (nargin < 7 | isempty(t))
|
||||
t = ones(m,1);
|
||||
end
|
||||
|
||||
if (nargin < 8 | isempty(mu))
|
||||
mu = zeros(m,1);
|
||||
end
|
||||
|
||||
t = full(t(:)); mu = full(mu(:));
|
||||
if (length(t) ~= m | length(mu) ~= m)
|
||||
error('Polyhedral input arguments are of incompatible sizes');
|
||||
end
|
||||
|
||||
l_p = -Big*ones(m,1);
|
||||
u_p = Big*ones(m,1);
|
||||
|
||||
idx = find(t > 0);
|
||||
if (length(idx) > 0)
|
||||
l_p(idx) = zeros(length(idx),1);
|
||||
end
|
||||
|
||||
idx = find(t < 0);
|
||||
if (length(idx) > 0)
|
||||
u_p(idx) = zeros(length(idx),1);
|
||||
end
|
||||
|
||||
mu = min(max(mu,l_p),u_p);
|
||||
else
|
||||
if (nargin >= 8 & ~isempty(mu))
|
||||
error('No polyhedral constraints -- multipliers set.');
|
||||
end
|
||||
|
||||
if (nargin >= 7 & ~isempty(t))
|
||||
error('No polyhedral constraints -- equation types set.');
|
||||
end
|
||||
end
|
||||
else
|
||||
A = [];
|
||||
end
|
||||
|
||||
% this is a fix, nnz may be bigger than this
|
||||
[f,J,domerr] = feval(cpfj,z+1e-5*ones(size(z))+1e-5*abs(z),1);
|
||||
|
||||
if (domerr > 0)
|
||||
[f,J,domerr] = feval(cpfj,z,1);
|
||||
end
|
||||
|
||||
if (domerr > 0)
|
||||
error([cpfj ' not defined at starting point']);
|
||||
end
|
||||
|
||||
if ~issparse(J)
|
||||
error([cpfj ' must return a sparse Jacobian']);
|
||||
end
|
||||
|
||||
nnzJ = nzmax(J);
|
||||
|
||||
row = n + m;
|
||||
ele = nnzJ + 2*nzmax(A);
|
||||
|
||||
init = [z; mu];
|
||||
low = [l; l_p];
|
||||
upp = [u; u_p];
|
||||
|
||||
if m > 0
|
||||
global mcp_vifunc;
|
||||
global mcp_viconn;
|
||||
global mcp_viconm;
|
||||
global mcp_viconA;
|
||||
global mcp_viconb;
|
||||
|
||||
mcp_vifunc = cpfj;
|
||||
mcp_viconn = n;
|
||||
mcp_viconm = m;
|
||||
mcp_viconA = A;
|
||||
mcp_viconb = b;
|
||||
|
||||
[status, ttime, f, J] = mcppath(row, ele, init, low, upp, 'mcp_vifunjac');
|
||||
else
|
||||
[status, ttime, f, J] = mcppath(row, ele, init, low, upp, cpfj);
|
||||
end
|
||||
|
||||
%if (status ~= 1)
|
||||
% status,
|
||||
% error('Path fails to solve problem');
|
||||
%end
|
||||
|
||||
mu = [];
|
||||
z = init;
|
||||
|
||||
if m > 0
|
||||
mu = init(n+1:n+m);
|
||||
z = init(1:n);
|
||||
|
||||
J = J(1:n,1:n);
|
||||
f = f(1:n) + A'*mu;
|
||||
end
|
||||
|
||||
return;
|
||||
Loading…
Reference in New Issue