Topic: algencan and ipopt via nls to nlsp convertor
Problem: to solve a non-linear system.
I was able to solve my system with 'scipy_fsolve', 'nssolve' and 'nlp:ralg'. Result is not too good, but seems to be true inside some limited area. Then I compiled algencan and ipopt+pyipopt from sources (and no errors were faced) and tried to use them, but I didn't succeed with either. The following errors occured:
ALGENCAN (line 112 in the source code):
-----------------------------------------------------
solver: algencan problem: unnamed
iter objFunVal
0 0.000e+00
Traceback (most recent call last):
File "old_solve_eqs.py", line 112, in <module>
r = p.solve('nlp:algencan') # still doesn't work
File "/usr/local/lib/python2.6/dist-packages/openopt-0.28-py2.6.egg/openopt/kernel/baseProblem.py", line 219, in solve
return runProbSolver(self, *args, **kwargs)
File "/usr/local/lib/python2.6/dist-packages/openopt-0.28-py2.6.egg/openopt/kernel/runProbSolver.py", line 185, in runProbSolver
R = converter(solverName, **solver_params)
File "/usr/local/lib/python2.6/dist-packages/openopt-0.28-py2.6.egg/openopt/kernel/NLSP.py", line 82, in nlsp2nlp
r = p.solve(solver, **solver_params)
File "/usr/local/lib/python2.6/dist-packages/openopt-0.28-py2.6.egg/openopt/kernel/baseProblem.py", line 219, in solve
return runProbSolver(self, *args, **kwargs)
File "/usr/local/lib/python2.6/dist-packages/openopt-0.28-py2.6.egg/openopt/kernel/runProbSolver.py", line 216, in runProbSolver
if not hasattr(p, 'xf') or all(p.xf==nan): p.xf = p.xk
AttributeError: NLP instance has no attribute 'xk'Ipopt (line 31 in the source code):
-----------------------------------------------------
solver: ipopt problem: unnamed
iter objFunVal log10(maxResidual)
0 0.000e+00 -100.00
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Common Public License (CPL).
For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************
NOTE: You are using Ipopt by default with the MUMPS linear solver.
Other linear solvers might be more efficient (see Ipopt documentation).
[Error] Ipopt faied in solving problem instance
Traceback (most recent call last):
File "old_solve_eqs.py", line 113, in <module>
r = p.solve('nlp:ipopt') # still doesn't work
File "/usr/local/lib/python2.6/dist-packages/openopt-0.28-py2.6.egg/openopt/kernel/baseProblem.py", line 219, in solve
return runProbSolver(self, *args, **kwargs)
File "/usr/local/lib/python2.6/dist-packages/openopt-0.28-py2.6.egg/openopt/kernel/runProbSolver.py", line 185, in runProbSolver
R = converter(solverName, **solver_params)
File "/usr/local/lib/python2.6/dist-packages/openopt-0.28-py2.6.egg/openopt/kernel/NLSP.py", line 82, in nlsp2nlp
r = p.solve(solver, **solver_params)
File "/usr/local/lib/python2.6/dist-packages/openopt-0.28-py2.6.egg/openopt/kernel/baseProblem.py", line 219, in solve
return runProbSolver(self, *args, **kwargs)
File "/usr/local/lib/python2.6/dist-packages/openopt-0.28-py2.6.egg/openopt/kernel/runProbSolver.py", line 190, in runProbSolver
solver(p)
File "/usr/local/lib/python2.6/dist-packages/openopt-0.28-py2.6.egg/openopt/solvers/CoinOr/ipopt_oo.py", line 175, in __solver__
x, zl, zu, obj = nlp.solve(p.x0)
TypeError: result of eval_f must be a floatProbably my code not too clean and I have to say that I have no complete understanding of FuncDesigner principles (BTW, where can I read the best manual?). But I suppose that if code works for some solvers, then it might be bugs in implementations of others?
Full source code is below.
Initial problem is inherited from physics (some crystalline nano-blah-blah-magnetism). Equations are quite complicated, but I can show them, if it is necessary.
Yury, Institute of Physics, Far-Eastern National University.
#coding: utf-8
import numpy as np
import scipy.spatial as sp
from numpy import abs, sign
""" If my own exceptions will be implemented"""
class LatticeError(Exception): pass
A = 1.001 # parameter of the lattice p — for NNS
n = 2
N = n**3
def GenerateCubicLattice (N):
""" Generate <cubic (N x N x N)> lattice """
lattice = np.arange(1,N**3+1).reshape(N,N,N)
return lattice
Lattice = GenerateCubicLattice(n)
"""Miller's indices: j --> (i,j,k)"""
def Miller(lattice, item):
if len(lattice.shape) == 1:
try:
return [lattice.tolist().index(item)+1]
except:
pass
else:
for i, subarray in enumerate(lattice):
try:
return [i+1] + Miller(subarray, item)
except:
pass
def GenerateLatticeKDtree ():
""" Generate kd-tree for NN searching """
data = [ Miller(Lattice,elem) for elem in np.arange(1,len(Lattice.flatten())+1) ]
tree = sp.KDTree(data)
return tree
Tree = GenerateLatticeKDtree()
def k_j (j):
""" $k_j$ is a set of closest neighbours of the node number $j$"""
k_j = Tree.query_ball_point(Miller(Lattice,j),A)
k_j = [n+1 for n in k_j]
k_j.remove(j)
k_j.sort()
return k_j
def L_j (j):
k_j_array = k_j(j)
L = np.array([ [ (-1)**int(n) for n in np.binary_repr(elem).zfill(len(k_j_array))] for elem in np.arange(2**len(k_j_array))])
w_j = L*k_j_array
return w_j
def Omega_j (j):
return [ reduce(lambda x,y: x+'*'+y ,[\
'a['+str(j-1)+']' if j>0 else '(1-a['+str(-j-1)+'])'\
for j in elem]) for elem in L_j(j) ]
def M_j (l,j):
return reduce(lambda x,y: x+('' if y[0]=='-' else '+')+y,\
[ ('-' if elem<0 else '')+'abs(2*a['+str(abs(elem)-1)+']-1)'\
for elem in L_j(j)[l] ])
def eq (j):
return reduce(lambda x,y: x+'+'+y, [ Omega_j(j)[l]+'*tanh(('+M_j(l,j)+')/t)' for l in range(len(Omega_j(j))) ])+'-(2*a['+str(j-1)+']-1)'
def equations():
return tuple(eq(i) for i in range(1,n**3+1))
#=========================
#print Lattice[0,0,1]
#print "Miller(k_2)=",Miller(Lattice,2)
#print "k_2 =",k_j(2)
#print L_j(2)
#print "Omega_2 =",Omega_j(2)
#print "M_j(3,2) =",M_j(3,2)
#print "eq(1)=", eq(1)
#=========================
from FuncDesigner import *
from openopt import NLSP
from numpy import arange, ones, empty
from numpy import tanh as npTanh
from numpy import cosh as npCosh
def tanh(inp):
if not isinstance(inp, oofun): return npTanh(inp)
return oofun(npTanh, input = inp, d = lambda x: Diag(1.0 / npCosh(x) ** 2))
#print equations()
a = oovar(size=N)
M = []
startPoint = {a: ones(N)}
for t in arange(0,2.1,0.3):
f = tuple(eval(eq(i)) for i in range(1,N+1))
p = NLSP(f, startPoint, maxIter=1e9, control=1e-15)
#r = p.solve('scipy_fsolve') # quite bad
#r = p.solve('nssolve') # best for now
#r = p.solve('nlp:ralg') # okay
#r = p.solve('nlp:algencan') # still doesn't work
r = p.solve('nlp:ipopt') # still doesn't work
## PRINTING
print "\n***** t = %f *****" %t#," startPoint = ",startPoint.values()[0]
print
m = 0
for i in range(n):
for j in range(n):
for k in range(n):
print('a%i%i%i: %f' %(i,j,k,a[m](r))),
m+=1
print
print
M.append((2*a(r)-1).sum()/N)
print "< m > =",M[-1]
startPoint = {a: a(r)}
print "< m > =", MP.S. Attachments doesn't work here, is it okay?