Here's an example of general constrained NL problem and solution by OO solver lincher
(x0-5)^2 + (x1-5)^2 + ... +(x149-5)^2 -> min
subjected to
# lb<= x <= ub:
x4 <= 4
8 <= x5 <= 15
# Ax <= b
x0+...+x149 >= 800
x10+x11 <= 9
# Aeq x = beq
x100+x101 = 8
# c(x) <= 0
2*x0^4-32 <= 0
x1^2+x2^2-8 <= 0
# h(x) = 0
1e6*(x[149]-1)**4 = 0
(x[148]-1.5)**4 = 0
----------------------------------------------------------------
from scikits.openopt import NLP
from numpy import cos, arange, ones, asarray, zeros
N = 150
p = NLP(lambda x: ((x-5)**2).sum(), 8*cos(arange(N)), iprint = 25, maxIter = 1e3)
#f(x) gradient (optional):
p.df = lambda x: 2*(x-5)
p.c = lambda x: [2* x[0] **4-32, x[1]**2+x[2]**2 - 8]
#c(x) gradients (optional):
nc = 2 # number of c(x)<=0 constraints
def DC(x):
r = zeros((p.n, nc))
r[0,0] = 2 * 4 * x[0]**3
r[1,1] = 2 * x[1]
r[2,1] = 2 * x[2]
return r
p.dc = DC
h1 = lambda x: 1e6*(x[-1]-1)**4
h2 = lambda x: (x[-2]-1.5)**4
p.h = [h1, h2]
#h(x) gradients (optional):
nh = len(p.h)
def DH(x):
r = zeros((p.n, nh))
r[-1,0] = 1e6*4*(x[-1]-1)**3
r[-2,1] = 4*(x[-2]-1.5)**3
return r
p.dh = DH
p.lb = -6*ones(p.n)
p.ub = 6*ones(p.n)
p.ub[4] = 4
p.lb[5], p.ub[5] = 8, 15
p.A = -ones((2, p.n))
p.A[1, 2:] = 0
p.A[1, 10:12] = 1
p.b = [-800, 9]
p.Aeq = zeros(p.n)
p.Aeq[100:102] = 1
p.beq = 8
r = p.solve('lincher')
--------------------------------------------
>>> starting solver lincher (BSD license) with problem unnamed
itn 0: Fk= 8596.39550577 maxResidual= 60585923.7208
itn 25 : Fk= 99.1913689183 maxResidual= 0.219982291475
N= 5534.59379932 alpha= 0.999266862564
itn 50 : Fk= 98.8504029574 maxResidual= 0.0025405684671
N= 111496.623829 alpha= 0.41746278006
itn 70 : Fk= 98.8329244399 maxResidual= 9.99877777019e-07
N= 115603.527987 alpha= 0.424851108558
solver lincher has finished solving the problem unnamed
istop: 4 (|| F[k] - F[k-1] || < funtol)
Solver: Time elapsed = 7.25 CPU Time Elapsed = 6.95
objFunValue: 98.8329244399 (feasible)
Friday, June 29, 2007
Wednesday, June 27, 2007
NLP solver "lincher": already works
lincher means "LINearization solver written in CHERkassy town, Ukraine"
currently it requires QP solver (the only one connected to OO is cvxopt one, also mosek)
Here's a small example, tomorrow I will have it made more exact:
from scikits.openopt import NLP
from numpy import cos, arange, ones, asarray
N = 15
p = NLP(lambda x: (x**2).sum()+((x-10)**2).sum(), cos(arange(N)))
p.c = lambda x: [2* x[0] **4-32, x[1]**2+x[2]**2 - 8]
h1 = lambda x: 1e6*(x[-1]-1)**4
h2 = lambda x: (x[-2]-1.5)**4
p.h = [h1, h2]
p.lb = -6*ones(p.n)
p.ub = 6*ones(p.n)
p.maxIter = 1e3
r = p.solve('lincher')
>>> starting solver lincher (BSD license) with problem unnamed
itn 0: Fk= 1485.60535217 maxResidual= 555356.768901
itn 25 : Fk= 861.006426698 maxResidual= 0.0815128971501
N= 196862.294712 alpha= 0.0140694378966
itn 50 : Fk= 860.925453409 maxResidual= 0.000180521169415
N= 197037.161793 alpha= 0.155520860502
itn 75 : Fk= 860.890756449 maxResidual= 1.74602521597e-06
N= 200312.271054 alpha= 0.0503775225937
itn 85 : Fk= 860.887158576 maxResidual= 9.89853580435e-07
N= 203978.915945 alpha= 0.0328485602443
solver lincher finished solving the problem unnamed
istop: 4 (|| F[k] - F[k-1] || < funtol)
Solver: Time elapsed = 1.49 CPU Time Elapsed = 1.42
objFunValue: 860.887158576 (feasible)
>>> r.xf
array([ 2.00000001, 2.00015789, 1.99984234, 5.00000007, 5.0000001 ,
5.0000001 , 4.99999931, 5.00000009, 5.00000007, 5.00000008,
5.00000002, 4.99999999, 4.99999994, 1.47340852, 0.99915843])
other changes:
r.msg field has been added (see line above:
istop: 4 (|| F[k] - F[k-1] || < funtol)
)
r contains iterFvals and iterXvals - values of x[k] and f[k] - k-th iter X and objFunc(x)
+some minor changes
currently it requires QP solver (the only one connected to OO is cvxopt one, also mosek)
Here's a small example, tomorrow I will have it made more exact:
from scikits.openopt import NLP
from numpy import cos, arange, ones, asarray
N = 15
p = NLP(lambda x: (x**2).sum()+((x-10)**2).sum(), cos(arange(N)))
p.c = lambda x: [2* x[0] **4-32, x[1]**2+x[2]**2 - 8]
h1 = lambda x: 1e6*(x[-1]-1)**4
h2 = lambda x: (x[-2]-1.5)**4
p.h = [h1, h2]
p.lb = -6*ones(p.n)
p.ub = 6*ones(p.n)
p.maxIter = 1e3
r = p.solve('lincher')
>>> starting solver lincher (BSD license) with problem unnamed
itn 0: Fk= 1485.60535217 maxResidual= 555356.768901
itn 25 : Fk= 861.006426698 maxResidual= 0.0815128971501
N= 196862.294712 alpha= 0.0140694378966
itn 50 : Fk= 860.925453409 maxResidual= 0.000180521169415
N= 197037.161793 alpha= 0.155520860502
itn 75 : Fk= 860.890756449 maxResidual= 1.74602521597e-06
N= 200312.271054 alpha= 0.0503775225937
itn 85 : Fk= 860.887158576 maxResidual= 9.89853580435e-07
N= 203978.915945 alpha= 0.0328485602443
solver lincher finished solving the problem unnamed
istop: 4 (|| F[k] - F[k-1] || < funtol)
Solver: Time elapsed = 1.49 CPU Time Elapsed = 1.42
objFunValue: 860.887158576 (feasible)
>>> r.xf
array([ 2.00000001, 2.00015789, 1.99984234, 5.00000007, 5.0000001 ,
5.0000001 , 4.99999931, 5.00000009, 5.00000007, 5.00000008,
5.00000002, 4.99999999, 4.99999994, 1.47340852, 0.99915843])
other changes:
r.msg field has been added (see line above:
istop: 4 (|| F[k] - F[k-1] || < funtol)
)
r contains iterFvals and iterXvals - values of x[k] and f[k] - k-th iter X and objFunc(x)
+some minor changes
Monday, June 25, 2007
QP class: ready!
QP: constructor for Quadratic Problem assignment
1/2 x' H x + f' x -> min
subjected to
A x <= b Aeq x = beq lb <= x <= ub Alternatively to A/Aeq you can use Awhole matrix as it's descriebed in LP documentation (or both A, Aeq, Awhole) examples of valid calls: p = QP(H, f, params)
p = QP(numpy.ones((3,3)), f=numpy.array([1,2,4]), params)
p = QP(f=range(8)+15, H = numpy.diag(numpy.ones(8)), params)
p = QP(H, f, A=A, Aeq=Aeq, b=b, beq=beq, lb=lb, ub=ub, params)
INPUT:
H: size n x n matrix, symmetric, positive-definite
f: size n x 1 vector
A: size m1 x n matrix, subjected to A * x <= b Aeq: size m2 x n matrix, subjected to Aeq * x = beq OUTPUT: OpenOpt QP class instance Solving of QPs is performed via r = p.solve(string_name_of_solver) Currently string_name_of_solver can be: cvxopt_qp (GPL) - requires CVXOPT (http://abel.ee.ucla.edu/cvxopt) (unfortunately, I can't suppress cvxopt text output, its solvers don't have appropriate params) Let's concider the problem x1^2 + 2x2^2 + 3x3^2 + 15x1 + 8x2 + 80x3 -> min (1)
subjected to
x1 + 2x2 + 3x3 <= 150 (2)
8x1 + 15x2 + 80x3 <= 800 (3)
x2 - x3 = 25 (4)
x1 <= 15 (5)
from numpy import diag, matrix
from scikits.openopt import QP
p = QP(diag([1,2,3]), [15,8,80], A = matrix('1 2 3; 8 15 80'), b = [150, 800], Aeq = [0, 1, -1], beq = 25, ub = [15,inf,inf])
# or p = QP(H=diag([1,2,3]), f=[15,8,80], A = matrix('1 2 3; 8 15 80'), b = [150, 800], Aeq = [0, 1, -1], beq = 25, ub = [15,inf,inf])
r = p.solve('cvxopt_qp')
f_opt, x_opt = r.ff, r.xf
# x_opt = array([-14.99999995, -2.59999996, -27.59999991])
# f_opt = -2453.7999916
1/2 x' H x + f' x -> min
subjected to
A x <= b Aeq x = beq lb <= x <= ub Alternatively to A/Aeq you can use Awhole matrix as it's descriebed in LP documentation (or both A, Aeq, Awhole) examples of valid calls: p = QP(H, f, params)
p = QP(numpy.ones((3,3)), f=numpy.array([1,2,4]), params)
p = QP(f=range(8)+15, H = numpy.diag(numpy.ones(8)), params)
p = QP(H, f, A=A, Aeq=Aeq, b=b, beq=beq, lb=lb, ub=ub, params)
INPUT:
H: size n x n matrix, symmetric, positive-definite
f: size n x 1 vector
A: size m1 x n matrix, subjected to A * x <= b Aeq: size m2 x n matrix, subjected to Aeq * x = beq OUTPUT: OpenOpt QP class instance Solving of QPs is performed via r = p.solve(string_name_of_solver) Currently string_name_of_solver can be: cvxopt_qp (GPL) - requires CVXOPT (http://abel.ee.ucla.edu/cvxopt) (unfortunately, I can't suppress cvxopt text output, its solvers don't have appropriate params) Let's concider the problem x1^2 + 2x2^2 + 3x3^2 + 15x1 + 8x2 + 80x3 -> min (1)
subjected to
x1 + 2x2 + 3x3 <= 150 (2)
8x1 + 15x2 + 80x3 <= 800 (3)
x2 - x3 = 25 (4)
x1 <= 15 (5)
from numpy import diag, matrix
from scikits.openopt import QP
p = QP(diag([1,2,3]), [15,8,80], A = matrix('1 2 3; 8 15 80'), b = [150, 800], Aeq = [0, 1, -1], beq = 25, ub = [15,inf,inf])
# or p = QP(H=diag([1,2,3]), f=[15,8,80], A = matrix('1 2 3; 8 15 80'), b = [150, 800], Aeq = [0, 1, -1], beq = 25, ub = [15,inf,inf])
r = p.solve('cvxopt_qp')
f_opt, x_opt = r.ff, r.xf
# x_opt = array([-14.99999995, -2.59999996, -27.59999991])
# f_opt = -2453.7999916
Friday, June 22, 2007
QP class is required
Most of NLP linearisation algorithms from ASAI are based on QP, not LP subproblems. I have tried some lp-based and I'm not fond of results - convergence isn't good enough.
So I have started to implement QP class (some code already have been added to svn; first of all I intend to connect cvxopt default qp solver). I had spent much time thinking what is the best letter for the QP matrix: MATLAB uses H, cvxopt - P, TOMLAB uses F, etc, etc.
I desided to chose H, because it seems the most pretty to me, however, the drawback of the latter - some people can associate the one with Hesse matrix (for example, TOMLAB uses the notation). In OO Hesse matrix will be passed as d2f (also I intend to implement d2c, d2h, d2L), for example:
p = NLP(..., c=myIneqConstr, dc=myDIneqConstr, d2c=myD2IneqConstr)
p.h = myEqConstr
p.dh = myDEqConstr
p.d2h = myD2EqConstr
(myIneqConstr etc are names of corresponding functions).
So I have started to implement QP class (some code already have been added to svn; first of all I intend to connect cvxopt default qp solver). I had spent much time thinking what is the best letter for the QP matrix: MATLAB uses H, cvxopt - P, TOMLAB uses F, etc, etc.
I desided to chose H, because it seems the most pretty to me, however, the drawback of the latter - some people can associate the one with Hesse matrix (for example, TOMLAB uses the notation). In OO Hesse matrix will be passed as d2f (also I intend to implement d2c, d2h, d2L), for example:
p = NLP(..., c=myIneqConstr, dc=myDIneqConstr, d2c=myD2IneqConstr)
p.h = myEqConstr
p.dh = myDEqConstr
p.d2h = myD2EqConstr
(myIneqConstr etc are names of corresponding functions).
Thursday, June 21, 2007
ALGENCAN: maybe, promising OO NLP solver?
As I have mentioned in suggestions about COIN-OR IPOPT and other solvers, adding several solvers with similar algorithms isn't a best idea, because they will suffer the same drawbacks. I guess if a user's vote would be organized "what do you want first of all from openopt", first of all they would chose "several good NLP solvers with as much types of constraints as possible, as well as 1st and 2nd derivatives", such as implementation of SQP, rSQP, IP, GRG etc (one from each class).
Now I'm busy with the OO native linearisation-based solver ('"lincher", today some more code have been added to svn repository), but maybe in future the University of Campinas and University of São Paulo (Brazilia) joint project ALGENCAN (GPL) solver should be added. The most important feature - it has Python bindings. I compiled the one rather easily in my Linux system (unfortunately, windows users can expect some problems) and executed a test from py-file example. However, some things are organised very inconvenient - ALGENCAN authors propose to solve a python-written optimization problem via modifying a py-file attached (toyprob.py), and some funcs there have rather non-standard interface.
Also, ALGENCAN misses linear constraints (Ax<=b, Aeq x = beq), there is only c(x)<=0, h(x)=0, their 1st derivatives (optionally) and lb-ub bounds. Of course, I can add linear constraints as non-lin (+ corresponding derivatives obtaining) automatically, like I do (primarily for ralg solver) in my openopt for matlab/octave package.
Now I'm busy with the OO native linearisation-based solver ('"lincher", today some more code have been added to svn repository), but maybe in future the University of Campinas and University of São Paulo (Brazilia) joint project ALGENCAN (GPL) solver should be added. The most important feature - it has Python bindings. I compiled the one rather easily in my Linux system (unfortunately, windows users can expect some problems) and executed a test from py-file example. However, some things are organised very inconvenient - ALGENCAN authors propose to solve a python-written optimization problem via modifying a py-file attached (toyprob.py), and some funcs there have rather non-standard interface.
Also, ALGENCAN misses linear constraints (Ax<=b, Aeq x = beq), there is only c(x)<=0, h(x)=0, their 1st derivatives (optionally) and lb-ub bounds. Of course, I can add linear constraints as non-lin (+ corresponding derivatives obtaining) automatically, like I do (primarily for ralg solver) in my openopt for matlab/octave package.
Tuesday, June 19, 2007
suggestions about COIN-OR IPOPT and other solvers
I have received a question from Alan Isaac about IPOPT, and I decided it's better to describe my opinion here.
>The IPopt solver (in C++) might be a nice addition to OpenOpt? https://projects.coin-or.org/Ipopt
I know about the IPOPT solver, as well as whole COIN-OR project, very well. BTW I heard there is very nice LP solver, that is used in the linearisation-based IPOPT solver (and it's one of the reasons why IPOPT works so good).
It is really very good solver, according to the Hans D. Mittelmann benchmark results:
ampl-nlp
mpec (with equality constraints)
The problems are:
1) I had seen only 1-2 python bindings to COIN-OR solvers. As far as I remember, it was PuLP (A Python Linear Programming modeler), and the way it was implemented is unsuitable for solving large-scale tasks (it demands its own GAMS-like syntax for describing LPs).
Other python binding (that I got from web search) was dated to 2001 or so, and hence is currently unsupported.
I'm enable to organize efficient python-C++ binding, especially to IPOPT, because it requires passing function handles, that is much more difficult than just simple matrices, as it is in LP or QP problems (btw COIN-OR QP solver isn't very good, see here). So, the problem of writing coin-or to python bridge is hard by itself, and keeping the code up-to-date is also very hard time-consuming task.
2) IPOPT, as well as whole COIN-OR project, has CPL license, that has copyleft statement.
3) I has already started writing a linearisation - based NLP solver, and it will have BSD license. Since it has similar to IPOPT algorithm (linearisation - based), it will have about the same tasks where it yields same results (either good or bad) as IPOPT.
Particularly, for linearisation-based solvers main problems are:
1) if some linearised constraints forms a linear-dependent system in a point x[k] (i.e. obtained in k-th iteration), then it's very hard for LP solvers to handle the situation, and it requires much more time (and cputime). The same problem arises if x_opt has same linear-dependent system of some linearised constraints.
2) Usually (if objfunc and NL constraints are not very costly) most of time/cputime in linearisation solvers is spend to solving of LPs (or QPs, if quadratic approach is used), hence good LP/QP solver is required, as I have mentioned above.
Since I've connected 3 LP solvers (they already had python bindings: glpk with cvxopt binding, cvxopt native lp solver and lp_solve with its own binding; also, non-free mosek lp solver with cvxopt binding is available, however, I didn't test the one properly). So I think number of free LP solvers is enough for now, and connecting COIN-OR LP solver shouldn't be primary task; the same to writing IPOPT or some other coin-or solvers binding.
>The IPopt solver (in C++) might be a nice addition to OpenOpt? https://projects.coin-or.org/Ipopt
I know about the IPOPT solver, as well as whole COIN-OR project, very well. BTW I heard there is very nice LP solver, that is used in the linearisation-based IPOPT solver (and it's one of the reasons why IPOPT works so good).
It is really very good solver, according to the Hans D. Mittelmann benchmark results:
ampl-nlp
mpec (with equality constraints)
The problems are:
1) I had seen only 1-2 python bindings to COIN-OR solvers. As far as I remember, it was PuLP (A Python Linear Programming modeler), and the way it was implemented is unsuitable for solving large-scale tasks (it demands its own GAMS-like syntax for describing LPs).
Other python binding (that I got from web search) was dated to 2001 or so, and hence is currently unsupported.
I'm enable to organize efficient python-C++ binding, especially to IPOPT, because it requires passing function handles, that is much more difficult than just simple matrices, as it is in LP or QP problems (btw COIN-OR QP solver isn't very good, see here). So, the problem of writing coin-or to python bridge is hard by itself, and keeping the code up-to-date is also very hard time-consuming task.
2) IPOPT, as well as whole COIN-OR project, has CPL license, that has copyleft statement.
3) I has already started writing a linearisation - based NLP solver, and it will have BSD license. Since it has similar to IPOPT algorithm (linearisation - based), it will have about the same tasks where it yields same results (either good or bad) as IPOPT.
Particularly, for linearisation-based solvers main problems are:
1) if some linearised constraints forms a linear-dependent system in a point x[k] (i.e. obtained in k-th iteration), then it's very hard for LP solvers to handle the situation, and it requires much more time (and cputime). The same problem arises if x_opt has same linear-dependent system of some linearised constraints.
2) Usually (if objfunc and NL constraints are not very costly) most of time/cputime in linearisation solvers is spend to solving of LPs (or QPs, if quadratic approach is used), hence good LP/QP solver is required, as I have mentioned above.
Since I've connected 3 LP solvers (they already had python bindings: glpk with cvxopt binding, cvxopt native lp solver and lp_solve with its own binding; also, non-free mosek lp solver with cvxopt binding is available, however, I didn't test the one properly). So I think number of free LP solvers is enough for now, and connecting COIN-OR LP solver shouldn't be primary task; the same to writing IPOPT or some other coin-or solvers binding.
an other one dayly report
Today I have elapsed some hours for searching an efficient spline toolkit, that allows usage of 2nd order (quadratic) splines. The was some compiling problems with module spline from scipy.sandbox, only some minutes ago I have installed it. Tomorrow I shall learn how it works (i.e. read documentation).
I have uploaded scanned book "Linearisation method" by Ukrainian academician Phenichniy: now it's available here (in russian, I have no English version)
I decided to name the new constrained solver "lincher" - LINearisation solver that had been written in CHERkasy town, Ukraine. Some solver code has been added to svn.
I have uploaded scanned book "Linearisation method" by Ukrainian academician Phenichniy: now it's available here (in russian, I have no English version)
I decided to name the new constrained solver "lincher" - LINearisation solver that had been written in CHERkasy town, Ukraine. Some solver code has been added to svn.
Monday, June 18, 2007
NL constrained problems
Some NLP constraints code has been added.
However, I still nedd some more time for to learn the NLP algorithms received from ASAI.
However, I still nedd some more time for to learn the NLP algorithms received from ASAI.
Thursday, June 14, 2007
QP: later, NLP: learning alg
I decided to make creating of QP class some time later: I should think twice of QP constructor interface. I mean even for LP there are 3 main approaches for assigning linear constraints:
1) A x <= b, Aeq x = beq (MATLAB, CVXOPT) (implemented in OO)
2) A x {less, equal, greater} b (glpk, lp_solve) (implemented in OO)
3) b_l <= A x <= b_u (TOMLAB) (maybe I will implement the one in OO later)
3rd way is more powerful - it allows to reduce A matrix in a factor up to 2 for some cases, but, on the other hand, creating of such matrices A, b_l, b_u is less convenient, and I guess very small amount of LP solvers is capable of taking benefites from the case.
So, a number of QP assignment ways is even greater.
I have to spend several days for investigating some versions of NLP linearisation algorithm that I have obtained from applied systems analysis institute (ASAI), by dr. Danilin. Unfortunately, I haven't English version so I can't discuss them in my blog or scipy mailing lists.
1) A x <= b, Aeq x = beq (MATLAB, CVXOPT) (implemented in OO)
2) A x {less, equal, greater} b (glpk, lp_solve) (implemented in OO)
3) b_l <= A x <= b_u (TOMLAB) (maybe I will implement the one in OO later)
3rd way is more powerful - it allows to reduce A matrix in a factor up to 2 for some cases, but, on the other hand, creating of such matrices A, b_l, b_u is less convenient, and I guess very small amount of LP solvers is capable of taking benefites from the case.
So, a number of QP assignment ways is even greater.
I have to spend several days for investigating some versions of NLP linearisation algorithm that I have obtained from applied systems analysis institute (ASAI), by dr. Danilin. Unfortunately, I haven't English version so I can't discuss them in my blog or scipy mailing lists.
Wednesday, June 13, 2007
MILP class created
1) LPSolve has been renamed to lpSolve (requires lp_solve installed)
2) MILP class created
Currently the only one MILP solver connected to OO is lpSolve (requires lp_solve installed)
see
from scikits.openopt import MILP
help(MILP)
2) MILP class created
Currently the only one MILP solver connected to OO is lpSolve (requires lp_solve installed)
see
from scikits.openopt import MILP
help(MILP)
Tuesday, June 12, 2007
Monday, June 11, 2007
latest svn changes
1) docs improved
see
from scikits.openopt import LP, NLP, NSP
help(LP)
help(NLP)
help(NSP)
2) some changes in linear constraints.
Next things to do:
create MILP, QP classes, connect lp_solve to MILP and CVXOPT QP solver to openopt QP class
see
from scikits.openopt import LP, NLP, NSP
help(LP)
help(NLP)
help(NSP)
2) some changes in linear constraints.
Next things to do:
create MILP, QP classes, connect lp_solve to MILP and CVXOPT QP solver to openopt QP class
Friday, June 8, 2007
OpenOpt LP solvers
some days ago I published a short benchmark of LP solvers available in OpenOpt.
Currently they are:
LPSolve (LGPL) - requires lpsolve + Python bindings installations (all mentioned is available in here)
cvxopt_lp (GPL) - requires CVXOPT
cvxopt_glpk(GPL2) - requires CVXOPT & glpk
cvxopt_mosek(non-free) - requires CVXOPT & mosek
However, I have some problems with mosek (even with trial license), it don't want to run in my AMD 64 + Linux, yielding 'inner mosek error'. Maybe, you will be more lucky:)
So, that report showed LPSolve to be much more better than other solvers.
Now I found the matter (and it's not my bug).
So, this is due to handling of lb-ub bounds: they are absent in CVXOPT (hence in glpk, mosek), and I used a routine that translated lb-ub bounds to Ax> nVars (however, other tasks can yield other results).
BTW, I think that LP solvers by themself should handle lines with single non-zero, but these ones involved don't. So, in future (but I guess not nearest) I intend to add the handling to the OpenOpt Kernel.
However, it has disadvantage: no correct duals will be returned, or it will require restart LP solver with unchanged LP problem from x_opt once again.
--------------------------------------------------------------------------------------
See whole LP bench results in my comment.
Currently they are:
LPSolve (LGPL) - requires lpsolve + Python bindings installations (all mentioned is available in here)
cvxopt_lp (GPL) - requires CVXOPT
cvxopt_glpk(GPL2) - requires CVXOPT & glpk
cvxopt_mosek(non-free) - requires CVXOPT & mosek
However, I have some problems with mosek (even with trial license), it don't want to run in my AMD 64 + Linux, yielding 'inner mosek error'. Maybe, you will be more lucky:)
So, that report showed LPSolve to be much more better than other solvers.
Now I found the matter (and it's not my bug).
So, this is due to handling of lb-ub bounds: they are absent in CVXOPT (hence in glpk, mosek), and I used a routine that translated lb-ub bounds to Ax
BTW, I think that LP solvers by themself should handle lines with single non-zero, but these ones involved don't. So, in future (but I guess not nearest) I intend to add the handling to the OpenOpt Kernel.
However, it has disadvantage: no correct duals will be returned, or it will require restart LP solver with unchanged LP problem from x_opt once again.
--------------------------------------------------------------------------------------
See whole LP bench results in my comment.
Thursday, June 7, 2007
hunting for bugs: history repeats
Almost whole day (today) I spent hunting for ralg() bug. It yields difference with OpenOpt for MATLAB (with my test case, where I noticed incorrect behaviour) only after 4th iteration. When I finally fix the one, I remembered that I had spent similar time for same bug in my OpenOpt for MATLAB (year ago) and several hours in OpenOpt for fortress (some months ago) (there are nothing but the ralg for now).
OpenOpt integrated to scikits svn
At last!...
I can't say all is perfect, but, at least, it works.
Some code were stolen from other scikits projects (pymat, pyaudiolab) setup files.
I noticed David Cournapeau uses "... = i for i in ...", but it is available only since Python 2.5, while scipy developers still think should they keep support for 2.3. I had "i for i in ..." in my code but after the message about 2.3 in mailing lists I had commented the lines for a while.
Some days have been spent for studing the svn and setuptools stuff.
One more thing that I dislike in scikits svn is repeating module name:
as you see, 'mlabwrap' is present for twice
also, __init__ can be omitted (since Python 2.5), and I'm not fond of creating and maintaining of additional files.
--------------------------------------------------------------------------------------
Ok, so now I intend to spend some days adding MILP and QP classes to OpenOpt.
I can't say all is perfect, but, at least, it works.
Some code were stolen from other scikits projects (pymat, pyaudiolab) setup files.
I noticed David Cournapeau uses "... = i for i in ...", but it is available only since Python 2.5, while scipy developers still think should they keep support for 2.3. I had "i for i in ..." in my code but after the message about 2.3 in mailing lists I had commented the lines for a while.
Some days have been spent for studing the svn and setuptools stuff.
One more thing that I dislike in scikits svn is repeating module name:
trunk(taken from scikits page)
mlabwrap/
setup.py # setuptools based
setup.cfg # optional, needed e.g. for development releases
README.txt
scikits/
__init__.py # __import__('pkg_resources').declare_namespace(__name__)
mlabwrap/
as you see, 'mlabwrap' is present for twice
also, __init__ can be omitted (since Python 2.5), and I'm not fond of creating and maintaining of additional files.
--------------------------------------------------------------------------------------
Ok, so now I intend to spend some days adding MILP and QP classes to OpenOpt.
Tuesday, June 5, 2007
OpenOpt for Python
Hi all,
this blog is related to OpenOpt module for scikits (see scipy.org) - free Python-based optimization alternative to TOMLAB, AMPL, GAMS etc.
I had some experience with OpenOpt for MATLAB/Octave before. You can take a look at reviews and download that one from here
Now Python version includes 2 non-smooth UC solvers: ralg (an implementation of Naum Z.Shor r-alg) and ShorEllipsoid (Naum Z.Shor modificated method of ellipsoids). Former is intended for medium-scale problems and is used very intensively in lots software solutions of our non-smooth optimization department, cybernetics institute, Ukrainian science academy; latter is very special solver for small-scale problems only (and knowing r0 required), I add the one only because it's very simple to implement.
There is basic graphic output implemented (provided you have matplotlib installed), see example.py
Also, connections to some LP solvers have been made, see exampleLP.py.
I shall connect openopt to scikits svn during 1-2 days (I need to decide how to split all my files into subdirectories very carefully and learn how to use svn - previously I used only cvs). Currently you can download the OpenOpt for Python from here
WBR, Dmitrey
Cherkassy, Ukraine
this blog is related to OpenOpt module for scikits (see scipy.org) - free Python-based optimization alternative to TOMLAB, AMPL, GAMS etc.
I had some experience with OpenOpt for MATLAB/Octave before. You can take a look at reviews and download that one from here
Now Python version includes 2 non-smooth UC solvers: ralg (an implementation of Naum Z.Shor r-alg) and ShorEllipsoid (Naum Z.Shor modificated method of ellipsoids). Former is intended for medium-scale problems and is used very intensively in lots software solutions of our non-smooth optimization department, cybernetics institute, Ukrainian science academy; latter is very special solver for small-scale problems only (and knowing r0 required), I add the one only because it's very simple to implement.
There is basic graphic output implemented (provided you have matplotlib installed), see example.py
Also, connections to some LP solvers have been made, see exampleLP.py.
I shall connect openopt to scikits svn during 1-2 days (I need to decide how to split all my files into subdirectories very carefully and learn how to use svn - previously I used only cvs). Currently you can download the OpenOpt for Python from here
WBR, Dmitrey
Cherkassy, Ukraine
Subscribe to:
Posts (Atom)
