Free Python optimization framework

Friday, June 29, 2007

constrained NLP example with gradients

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)

1 comment:

Dmitrey said...

this example is obsolete, see this one: http://projects.scipy.org/scipy/scikits/browser/trunk/openopt/scikits/openopt/examples/nlp_1.py