Optimizing Fortran function in pyomo

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 -

How to successfully solve this uncostrained optimization problem using pyomo?
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?
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 -

Pyomo - 5.5.1 (VOTD) (CPython 2.7.12 on Linux 4.4.0-127-generic)
ipopt - Ipopt 3.12.8
f2py - installed along with numpy 1.16.2

Thanks for your help!

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!