Optimizing Fortran function in pyomo

Question

I wish to optimize a Fortran function using Pyomo. Both the objective function and the constraints are are written in Fortran. Based on the answer given here, we can use ExternalFunction expression object. But I am not able to get the results even for the simplest function. Given below is a reproducible example which consists of the Fortran function, the python (Python 2.7.12) script, the commands executed for optimization and the error.

Fortran function file (funcs.f) -

  SUBROUTINE OBJ1(ARG,OBJ)
  DOUBLE PRECISION, INTENT(IN) :: ARG(2)
  DOUBLE PRECISION, INTENT(OUT) :: OBJ
  OBJ = ARG(1)+ARG(2)
  RETURN
  END SUBROUTINE

Python script (pytest.py) -

import funcs
from pyomo.environ import *
from pyomo.opt import *
from pyomo.core import *

m = ConcreteModel()
m.a = Var(RangeSet(1,2),within=NonNegativeReals,bounds=(0,10))
m.f = ExternalFunction(library='funcs.so',function='OBJ1')
expr = m.f(m.a)
m.obj = Objective(expr=expr,sense=minimize)

opt = SolverFactory('ipopt')
results = opt.solve(m,tee=True)
print(results)

Commands executed at the terminal -

>> f2py -c -m funcs funcs.f
>> python pytest.py

Error -

File "/usr/local/lib/python2.7/dist-packages/pyomo/core/base/external.py", line 160, in load_library
FUNCADD(('funcadd_ASL', self._so))(byref(AE))
AttributeError: /home/utkarsh/Desktop/python/modules/blackboxOptimization/funcs.so: undefined symbol: funcadd_ASL

I have given only small portion of the error which I thought was relevant.

Given this, I have a the following questions -

1. How to successfully solve this uncostrained optimization problem using pyomo?

2. For my complete project, I will have to give constraints in Fortran itself. The constraint subroutines will return a real number which will be bounded using pyomo. How to model these type of constraints?

3. I am assuming that Pyomo takes this blackbox as non-linear optimization. Hence, I am using ipopt solver. Is this assumption correct?

The versions of packages -

1. Pyomo - 5.5.1 (VOTD) (CPython 2.7.12 on Linux 4.4.0-127-generic)

2. ipopt - Ipopt 3.12.8

3. f2py - installed along with numpy 1.16.2

Thanks for your help!

Answer 1

If you are not bound to Pyomo you could use the excellent Pygmo package which offers solvers for different kinds of problems including blackbox solvers.

Here's a small example on how to use it on a continuous constrained single objective test problem:

import pygmo as pg
import pandas as pd


class Rosenbrock():
    """Rosenbrock function constrained to a disk.

    See: https://en.wikipedia.org/wiki/Test_functions_for_optimization
    """

    def fitness(self, x):
        """Evaluate fitness.

        Instead of the Rosenbrock function you could call your Fortran
        code here e.g. by using F2PY: https://www.numfys.net/howto/F2PY/
        """
        obj = (1-x[0])**2+100*(x[1]-x[0]**2)**2
        ineq = x[0]**2+x[1]**2-2
        return [obj, ineq]

    def get_bounds(self):
        """Return boundaries."""
        return ([-1.5]*2, [1.5]*2)

    def get_nic(self):
        """Determine number of inequalities."""
        return 1


# set up and solve problem
pro = pg.problem(Rosenbrock())
pop = pg.population(pro, size=200)
# see: https://github.com/esa/pagmo2/blob/master/include/pagmo/algorithms/
algo = pg.algorithm(pg.ihs(gen=10000))
algo.set_verbosity(100)
pop = algo.evolve(pop)

# extract solutions
fits = pd.DataFrame(pop.get_f())
vectors = pd.DataFrame(pop.get_x())
best_idx = pop.best_idx()
best_vector = vectors.loc[best_idx].to_frame().T
best_fitness = fits.loc[best_idx].to_frame().T

print(best_vector)
print(best_fitness)

You would then just have to deal with "interfacing" your Fortran code within the fitness function.

Hope this helps!