Commit abbc0209 authored by Laurent Belcour's avatar Laurent Belcour

Working qpas solver on 1d and 2d generated data

parent 74ae3659
+----------------------------------------------+
| Written by Adrian Wills, |
| School of Elec. Eng. & Comp. Sci. |
| University of Newcastle, |
| Callaghan, NSW, 2308, AUSTRALIA |
| |
| Last Revised 18 January 2010. |
| |
| Copyright (C) Adrian Wills. |
+----------------------------------------------+
The current version of this software is free of charge and
openly distributed, BUT PLEASE NOTE:
This software must be referenced when used in a published work.
This software may not be re-distributed as a part of any product.
This software 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.
IF IT FAILS TO WORK, IT'S YOUR LOSS AND YOUR PROBLEM.
------------------------------------------------------------------------ +----------------------------------------------+ | Written by Adrian Wills, | | School of Elec. Eng. & Comp. Sci. | | University of Newcastle, | | Callaghan, NSW, 2308, AUSTRALIA | | | | Last Revised 11 August 2009. | | | | Copyright (C) Adrian Wills. | +----------------------------------------------+ The current version of this software is free of charge and openly distributed, BUT PLEASE NOTE: This software must be referenced when used in a published work. This software may not be re-distributed as a part of a product. This software 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. IF IT FAILS TO WORK, IT'S YOUR LOSS AND YOUR PROBLEM. ---------------------------------------------------------------------------
To get the software working, proceed as follows: 1. Once unzipped (you have probably done this bit already :o)) the archive should contain the following files: ./qpas.m ./qpas.mexmaci ./qpas.mexw32 ./qpas.mexglx ./qpas.mexa64 ./qpip.m ./qpip.mexmaci ./qpip.mexw32 ./qpip.mexglx ./qpip.mexa64 ./qps_as.m ./qps_as.mexmaci ./qps_as.mexw32 ./qps_as.mexglx ./qps_as.mexa64 ./qps_ip.mexmaci ./qps_ip.mexw32 ./qps_ip.mexglx ./qps_ip.mexa64 ./qps_mq.m ./qps_mq.mexmaci ./qps_mq.mexw32 ./qps_mq.mexglx ./qps_mq.mexa64 ./qpgen.m ./test.m ./README ./LICENSE 2. Move these files to a directory of your choice. 3. Start Matlab and change to the directory you used in step 2. 4. Try running the test.m script. That's it! Please email any questions/comments/criticisms to: Adrian.Wills@newcastle.edu.au
\ No newline at end of file
% This dual active-set algorithm attempts to solve
% quadratic programming problems in the following form:
%
% x* = arg min 0.5x'Hx + f'x
% x
%
% s.t. Ax = b, linear equality constraints,
% Lx <= k, general linear inequality constraints,
% l <= x <= u, bound constraints.
%
% In Matlab a call to this function would look like:
%
% [x,err,lm] = qp(H,f,L,k,A,b,l,u,display);
%
% where the inputs are:
%
% H: An (n x n) positive semi-definite symmetric matrix
%
% f: A n element column vector
%
% (L,k): General linear inequality constraints
%
% (A,b): General linear equality constraints
%
% l: Element-wise lower bound constraints
%
% u: Element-wise upper bound constraints
%
% display: If display>0 then iteration information is displayed
%
%
% and the outputs are:
%
% x: the optimal solution (if obtained)
%
% err: error number, if err=0, then x is optimal
%
% lm: structure of Lagrange multipliers
%
% lm.equality: Lagrange multipliers for equality constraints
%
% lm.inequality: Lagrange multipliers for general inequality constraints
%
% lm.lowerbound: Lagrange multipliers for lower bound constraints
%
% lm.upperbound: Lagrange multipliers for upper bound constraints
%
% If (A,b) and/or (L,k) and/or l and/or u, are not being used, then
% set the respective entry to [] (i.e. the empty matrix).
%
% E.g.1 suppose that only L and k are being used then a call would look
% like
%
% x = qp(H,f,L,k);
%
% E.g.2 suppose that (A,b) and l are required, then a call would look
% like
%
% x = qp(H,f,[],[],A,b,l);
%
% +----------------------------------------------+
% | Written by Adrian Wills, |
% | School of Elec. Eng. & Comp. Sci. |
% | University of Newcastle, |
% | Callaghan, NSW, 2308, AUSTRALIA |
% | |
% | Last Revised 25 May 2007. |
% | |
% | Copyright (C) Adrian Wills. |
% +----------------------------------------------+
%
% The current version of this software is free of charge and
% openly distributed, BUT PLEASE NOTE:
%
% This software must be referenced when used in a published work.
%
% This software may not be re-distributed as a part of a commercial product.
% If you distribute it in a non-commercial products, please contact me first,
% to make sure you ship the most recent version.
%
% This software 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.
%
% IF IT FAILS TO WORK, IT'S YOUR LOSS AND YOUR PROBLEM.
\ No newline at end of file
function [H,f,A,b,L,k,l,u] = qpgen(nx,ne,nc,nb,rank,def,z);
% This function generates a random quadratic program in the form:
%
% min 0.5x'Hx + f'x
% x
% s.t. Ax = b
% Lx <= k
% l <= x <= u.
%
% Usage is as follows:
%
% [H,f,A,b,L,k,l,u] = qpgen(nx,ne,nc,nb,rank,def,z);
%
% where nx is the dimension of x, ne is the number of equality constraints,
% nc is the number of general inequality constraints and nb is either equal
% to nx or is zero.
%
% Rank referes to the rank of H and def is 0 for positive semi-definite,
% and 1 for negative semi-definite. If z is 1 then H will contain zeros.
%
% E.g. copy the following to the Matlab command line to generate a convex
% quadratic programming problem with linear equality, linear inequality
% and simple bound constraints.
%
% n=100; [H,f,A,b,L,k,l,u]=qpgen(n,round(0.5*n),round(1.5*n),n,n,0,0);
%
% +----------------------------------------------+
% | Written by Adrian Wills, |
% | School of Elec. Eng. & Comp. Sci. |
% | University of Newcastle, |
% | Callaghan, NSW, 2308, AUSTRALIA |
% | |
% | Last Revised 25 May 2007. |
% | |
% | Copyright (C) Adrian Wills. |
% +----------------------------------------------+
%
% The current version of this software is free of charge and
% openly distributed, BUT PLEASE NOTE:
%
% This software must be referenced when used in a published work.
%
% This software may not be re-distributed as a part of a commercial product.
% If you distribute it in a non-commercial products, please contact me first,
% to make sure you ship the most recent version.
%
% This software 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.
%
% IF IT FAILS TO WORK, IT'S YOUR LOSS AND YOUR PROBLEM.
nx = max(nx,0);
ne = max(ne,0);
nc = max(nc,0);
nb = max(nx,0);
H = []; f = [];
L = []; k = [];
A = []; b = [];
l = []; u = [];
rank = max(rank,0);
rank = min(rank,nx);
if (def > 1 & def < 0), def = 0; end
if z ==1,
H = [randn(rank,rank);zeros(nx-rank,rank)];
else
H = randn(nx,rank);
end
if def == 0,
H = H*H';
else
H = (H+H')/2;
end
f = 100*randn(nx,1);
if nc,
L = randn(nc,nx);
k = rand(nc,1);
end
if ne,
A = randn(ne,nx);
b = rand(ne,1);
end
if nb,
l = -rand(nb,1);
u = rand(nb,1);
end
% This primal-dual predictor-corrector algorithm attempts to solve
% quadratic programming problems in the following form:
%
% x* = arg min 0.5x'Hx + f'x
% x
%
% s.t. Ax = b, linear equality constraints,
% Lx <= k, general linear inequality constraints,
% l <= x <= u, bound constraints.
%
% In Matlab a call to this function would look like:
%
% [x,err,lm] = qpip(H,f,L,k,A,b,l,u,display,mu,method);
%
% where the inputs are:
%
% H: An (n x n) positive semi-definite symmetric matrix
%
% f: A n element column vector
%
% (L,k): General linear inequality constraints
%
% (A,b): General linear equality constraints
%
% l: Element-wise lower bound constraints
%
% u: Element-wise upper bound constraints
%
% display: If display>0 then iteration information is displayed
%
% mu: Desired complementarity gap target (point on central-path)
% [Default is mu=0.0]
%
% method: If method=1 (default), then a faster but less accurate linear
% solve step is used. Conversely, if method=0 then a slower but
% more accurate linear solve step is used.
%
%
% and the outputs are:
%
% x: the optimal solution (if obtained)
%
% err: error number, if err=0, then x is optimal
%
% lm: structure of Lagrange multipliers
%
% lm.equality: Lagrange multipliers for equality constraints
%
% lm.inequality: Lagrange multipliers for general inequality constraints
%
% lm.lowerbound: Lagrange multipliers for lower bound constraints
%
% lm.upperbound: Lagrange multipliers for upper bound constraints
%
% If (A,b) and/or (L,k) and/or l and/or u, are not being used, then
% set the respective entry to [] (i.e. the empty matrix).
%
% E.g.1 suppose that only L and k are being used then a call would look
% like
%
% x = qpip(H,f,L,k);
%
% E.g.2 suppose that (A,b) and l are required, then a call would look
% like
%
% x = qpip(H,f,[],[],A,b,l);
%
% +----------------------------------------------+
% | Written by Adrian Wills, |
% | School of Elec. Eng. & Comp. Sci. |
% | University of Newcastle, |
% | Callaghan, NSW, 2308, AUSTRALIA |
% | |
% | Last Revised 25 May 2007. |
% | |
% | Copyright (C) Adrian Wills. |
% +----------------------------------------------+
%
% The current version of this software is free of charge and
% openly distributed, BUT PLEASE NOTE:
%
% This software must be referenced when used in a published work.
%
% This software may not be re-distributed as a part of a commercial product.
% If you distribute it in a non-commercial products, please contact me first,
% to make sure you ship the most recent version.
%
% This software 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.
%
% IF IT FAILS TO WORK, IT'S YOUR LOSS AND YOUR PROBLEM.
\ No newline at end of file
% This dual active-set algorithm attempts to solve
% quadratic programming problems in the following form:
%
% x* = arg min 0.5x'Hx + f'x
% x
%
% s.t. l <= x <= u, bound constraints.
%
% In Matlab a call to this function would look like:
%
% [x,err,lm] = qps_as(H,f,l,u,display);
%
% where the inputs are:
%
% H: An (n x n) positive semi-definite symmetric matrix
%
% f: A n element column vector
%
% l: Element-wise lower bound constraints
%
% u: Element-wise upper bound constraints
%
% display: If display>0 then iteration information is displayed
%
%
% and the outputs are:
%
% x: the optimal solution (if obtained)
%
% err: error number, if err=0, then x is optimal
%
% lm: structure of Lagrange multipliers
%
% lm.lowerbound: Lagrange multipliers for lower bound constraints
%
% lm.upperbound: Lagrange multipliers for upper bound constraints
%
% NOTE: Both l and u must be provided, however, +/-inf entries are
% acceptable so that "free" variables are catered for.
%
% +----------------------------------------------+
% | Written by Adrian Wills, |
% | School of Elec. Eng. & Comp. Sci. |
% | University of Newcastle, |
% | Callaghan, NSW, 2308, AUSTRALIA |
% | |
% | Last Revised 25 May 2007. |
% | |
% | Copyright (C) Adrian Wills. |
% +----------------------------------------------+
%
% The current version of this software is free of charge and
% openly distributed, BUT PLEASE NOTE:
%
% This software must be referenced when used in a published work.
%
% This software may not be re-distributed as a part of a commercial product.
% If you distribute it in a non-commercial products, please contact me first,
% to make sure you ship the most recent version.
%
% This software 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.
%
% IF IT FAILS TO WORK, IT'S YOUR LOSS AND YOUR PROBLEM.
\ No newline at end of file
% This primal-dual predictor-corrector algorithm attempts to solve
% quadratic programming problems in the following form:
%
% x* = arg min 0.5x'Hx + f'x
% x
%
% s.t. l <= x <= u, bound constraints.
%
% In Matlab a call to this function would look like:
%
% [x,err,lm] = qps_ip(H,f,l,u,display);
%
% where the inputs are:
%
% H: An (n x n) positive semi-definite symmetric matrix
%
% f: A n element column vector
%
% l: Element-wise lower bound constraints
%
% u: Element-wise upper bound constraints
%
% display: If display>0 then iteration information is displayed
%
%
% and the outputs are:
%
% x: the optimal solution (if obtained)
%
% err: error number, if err=0, then x is optimal
%
% lm: structure of Lagrange multipliers
%
% lm.lowerbound: Lagrange multipliers for lower bound constraints
%
% lm.upperbound: Lagrange multipliers for upper bound constraints
%
% NOTE: Both l and u must be provided, however, +/-inf entries are
% acceptable so that "free" variables are catered for.
%
% +----------------------------------------------+
% | Written by Adrian Wills, |
% | School of Elec. Eng. & Comp. Sci. |
% | University of Newcastle, |
% | Callaghan, NSW, 2308, AUSTRALIA |
% | |
% | Last Revised 25 May 2007. |
% | |
% | Copyright (C) Adrian Wills. |
% +----------------------------------------------+
%
% The current version of this software is free of charge and
% openly distributed, BUT PLEASE NOTE:
%
% This software must be referenced when used in a published work.
%
% This software may not be re-distributed as a part of a commercial product.
% If you distribute it in a non-commercial products, please contact me first,
% to make sure you ship the most recent version.
%
% This software 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.
%
% IF IT FAILS TO WORK, IT'S YOUR LOSS AND YOUR PROBLEM.
\ No newline at end of file
% This branch-and-bound type algorithm attempts to solve
% quadratic programming problems in the following form:
%
% x* = arg min 0.5x'Hx + f'x
% x
%
% s.t. l <= x <= u, bound constraints.
%
% In Matlab a call to this function would look like:
%
% [x,err,lm] = qps_mq(H,f,l,u);
%
% where the inputs are:
%
% H: An (n x n) positive semi-definite symmetric matrix
%
% f: A n element column vector
%
% l: Element-wise lower bound constraints
%
% u: Element-wise upper bound constraints
%
% and the outputs are:
%
% x: the optimal solution (if obtained)
%
% err: error number, if err=0, then x is optimal
%
% lm: structure of Lagrange multipliers
%
% lm.lowerbound: Lagrange multipliers for lower bound constraints
%
% lm.upperbound: Lagrange multipliers for upper bound constraints
%
% NOTE: Both l and u must be provided, however, +/-inf entries are
% acceptable so that "free" variables are catered for.
%
% +----------------------------------------------+
% | Written by Adrian Wills, |
% | School of Elec. Eng. & Comp. Sci. |
% | University of Newcastle, |
% | Callaghan, NSW, 2308, AUSTRALIA |
% | |
% | Last Revised 25 May 2007. |
% | |
% | Copyright (C) Adrian Wills. |
% +----------------------------------------------+
%
% The current version of this software is free of charge and
% openly distributed, BUT PLEASE NOTE:
%
% This software must be referenced when used in a published work.
%
% This software may not be re-distributed as a part of a commercial product.
% If you distribute it in a non-commercial products, please contact me first,
% to make sure you ship the most recent version.
%
% This software 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.
%
% IF IT FAILS TO WORK, IT'S YOUR LOSS AND YOUR PROBLEM.
\ No newline at end of file
% Run this file to test routines
%
% +----------------------------------------------+
% | Written by Adrian Wills, |
% | School of Elec. Eng. & Comp. Sci. |
% | University of Newcastle, |
% | Callaghan, NSW, 2308, AUSTRALIA |
% | |
% | Last Revised 25 May 2007. |
% | |
% | Copyright (C) Adrian Wills. |
% +----------------------------------------------+
%
% The current version of this software is free of charge and
% openly distributed, BUT PLEASE NOTE:
%
% This software must be referenced when used in a published work.
%
% This software may not be re-distributed as a part of a commercial product.
% If you distribute it in a non-commercial products, please contact me first,
% to make sure you ship the most recent version.
%
% This software 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.
%
% IF IT FAILS TO WORK, IT'S YOUR LOSS AND YOUR PROBLEM.
% Generate a strictly convex quadratic program (random)
n=40; [H,f,A,b,L,k,l,u]=qpgen(n,round(0.2*n),3*n,n,n,0,0);
% clear all
%
% load test_data
dsp=1;
% Test qp
disp(sprintf('\n'));
disp('Testing qpas routine');
try,
tic;xqpas=qpas(H,f,L,k,A,b,l,u,dsp);toc
disp('Test OK');
catch,
disp('There was an error');
end
% Test qpip
disp(sprintf('\n'));
disp('Testing qpip routine');
try,
tic;xqpip=qpip(H,f,L,k,A,b,l,u,dsp); toc
disp('Test OK');
catch,
disp('There was an error');
end
% Test qps_as
disp(sprintf('\n'));
disp('Testing qps_as routine');
try,
tic;xqps_as=qps_as(H,f,l,u,dsp);toc
disp('Test OK');
catch,
disp('There was an error');
end
% Test qps_ip
disp(sprintf('\n'));
disp('Testing qps_ip routine');
try,
tic;xqps_ip=qps_ip(H,f,l,u,dsp);toc
disp('Test OK');
catch,
disp('There was an error');
end
% Test qps_mq
disp(sprintf('\n'));
disp('Testing qps_mq routine');
try,
tic;xqps_mq=qps_mq(H,f,l,u);toc
disp('Test OK');
catch,
disp('There was an error');
end
disp(sprintf('\n'));
\ No newline at end of file
......@@ -271,13 +271,22 @@ bool rational_fitter_matlab::fit_data(const rational_data* d, int np, int nq, in
engEvalString(ep, "display(b)");
std::cout << output << std::endl ;
#endif
#ifdef USE_MATLAB
engEvalString(ep, "[x, fval, flag] = quadprog(H,f,A,b);");
#else
engEvalString(ep, "cd matlab;");
engEvalString(ep, "[x, err] = qpas(H,f,A,b);");
engEvalString(ep, "flag = err == 0.0;");
engEvalString(ep, "cd ..;");
#endif
#ifdef DEBUG
std::cout << output << std::endl ;
engEvalString(ep, "display(x)");
std::cout << output << std::endl ;
engEvalString(ep, "display(flag)");
std::cout << output << std::endl ;
#endif
mxDestroyArray(H);
mxDestroyArray(f);
mxDestroyArray(A);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment