Couenne solver via AMPL mod export

demo.py

from casadi import SX, vertcat, inf

# Declare variables
x = SX.sym("x")
y = SX.sym("y")
z = SX.sym("z")

# Formulate the NLP
f = x**2 + 100*z**2
g = z + (1-x)**2 - y
nlp = {'x':vertcat(x,y,z), 'f':f, 'g':g}

# Solve the Rosenbrock problem
x0  = [2.5,3.0,0.75]
lbg = [0]
ubg = [0]
ubx = [inf]*3
lbx = [inf]*3

integer = [False, False, False]

# To use locally through ampl, make sure you have the ampl executable folder on your PATH
#  The solver executable may also accept a generated .nl file
#
# To use on NEOS, submit generated .mod and .dat files ( https://neos-server.org/neos/solvers/index.html )
# For an overview of solvers, see http://plato.asu.edu/sub/pns.html

def solve_with_ampl_solver(name,solver_name,nlp,x0=None,p=None,lbg=None,ubg=None,lbx=None,ubx=None,integer=None):
  import casadi
  nlp = casadi.Function('nlp',nlp,["x","p"],["f","g"])
  nlp = nlp.expand("nlp", {"live_variables": False})

  nx = nlp.numel_in(0)
  np = nlp.numel_in(1)
  ng = nlp.numel_out(1)
  if x0 is None:
    x0 = [0]*nx
  if lbg is None:
    lbg = [-casadi.inf]*ng
  if ubg is None:
    ubg = [casadi.inf]*ng
  if lbx is None:
    lbx = [-casadi.inf]*nx
  if ubx is None:
    ubx = [casadi.inf]*nx
  if integer is None:
    integer = [False]*nx
  if p is None:
    p = casadi.DM(0, 1)




  instr = nlp.instructions_sx()

  nx = nlp.numel_in(0)
  np = nlp.numel_in(1)
  ng = nlp.numel_out(1)

  print(lbg.shape, ubg.shape, ng)

  with open(name+".mod","w") as mod:

    print("param p {1..%d};" % np,file=mod)
    print("param x0 {1..%d};" % nx,file=mod)
    print("param lbx {1..%d};" % nx,file=mod)
    print("param ubx {1..%d};" % nx,file=mod)
    print("param lbg {1..%d};" % ng,file=mod)
    print("param ubg {1..%d};" % ng,file=mod)
    for i in range(nx):
      if integer[i]:
        print("var x{ip} integer >= lbx[{ip}], <= ubx[{ip}] := x0[{ip}];".format(ip=i+1),file=mod)
      else:
        print("var x{ip} >= lbx[{ip}], <= ubx[{ip}] := x0[{ip}];".format(ip=i+1),file=mod)

    def helper(i):
      return "h" + str(i)


    tainted = {}

    # Loop over the algorithm
    for k in range(nlp.n_instructions()):

      # Get the atomic operation
      op = nlp.instruction_id(k)
      o = nlp.instruction_output(k)
      i = nlp.instruction_input(k)

      if(op==casadi.OP_CONST):
        print("param " + helper(o[0]) + " = %.18f" % nlp.instruction_constant(k)+ ";",file=mod)
      else:
        if op==casadi.OP_INPUT:
          if i[0]==0:
            print("var " + helper(o[0]) + " = x" + str(i[1]+1) +  ";",file=mod)
            tainted[o[0]] = True
          else:
            print("param " + helper(o[0]) + " = p[" + str(i[1]+1) +  "];",file=mod)
        elif op==casadi.OP_OUTPUT:
          print(output + " " + nlp.name_out(o[0]) + str(o[1]+1) +  " = " + helper(i[0]) + ";",file=mod)
        else:
          t = False
          for e in i:
            t = t or e in tainted
          output = "var" if t else "param"
          if t:
            tainted[o[0]] = True
          if op==casadi.OP_SQ:
            print(output + " " + helper(o[0]) + " = " + helper(i[0]) + "^2;",file=mod)
          elif op==casadi.OP_TWICE:
            print(output + " " + helper(o[0]) + " = " + helper(i[0]) + "*2;",file=mod)
          else:
            disp_in = [helper(a) for a in i]
            print(output + " " + helper(o[0]) + " = " + casadi.print_operator(instr[k],disp_in) + ";",file=mod)


    print("minimize f:  f1;",file=mod)

    for i in range(ng):
      print("subject to Constraints{ip}: \n  lbg[{ip}] <= g{ip} <= ubg[{ip}];".format(i=i,ip=i+1) ,file=mod)

  def out(a):
    if a == casadi.inf: return "Infinity"
    if -a == casadi.inf: return "-Infinity"
    return "%.18f" % a

  with open(name+".dat","w") as dat:
    print("data;", file=dat)
    print("param: x0 lbx ubx := ", file=dat)
    for i in range(nx):
      print("{ip} {x0} {lbx} {ubx}".format(ip=i+1,x0=out(x0[i]),lbx=out(lbx[i]),ubx=out(ubx[i])), file=dat)
    print(";", file=dat)
    print("param: lbg ubg := ", file=dat)
    for i in range(ng):
      print("{ip} {lbg} {ubg}".format(ip=i+1,lbg=out(lbg[i]),ubg=out(ubg[i])), file=dat)
    print(";", file=dat)

    print("param: p := ", file=dat)
    for i in range(np):
      print("{ip} {p}".format(ip=i+1,p=out(p[i])), file=dat)
    print(";", file=dat)

  with open(name+".run","w") as run:
    print("model {name}.mod".format(name=name),file=run)
    print("data {name}.dat".format(name=name),file=run)
    print("option solver "+solver_name+";",file=run)
    print("option show_stats 1;",file=run)
    print("option eexit -100000;",file=run)

    print("solve;",file=run)

    for i in range(nx):
      print("display x{ip} > {name}.out;".format(name=name,ip=i+1),file=run)

    print("close {name}.out;".format(name=name),file=run)

  import subprocess
  subprocess.Popen(["ampl","-og"+name,name+".mod",name+".dat"]).wait()
  subprocess.Popen(["ampl",name+".run"]).wait()

  xsol = [0]*nx
  with open(name+".out","r") as out:
    for l in out:
      if l.startswith("x"):
        var, val = l.split("=")
        xsol[int(var[1:])-1] = float(val)

  return xsol

xsol = solve_ampl("rosenbrock", nlp,x0=x0,lbx=lbx,ubx=ubx,lbg=lbg,ubg=ubg,integer=integer)

print("xsol", xsol)