Help with setting up an optimization problem with nested functions in pyomo

[Harpreet Singh]

Hi. I have been looking for an optimization package for a month now after I figured scipy.optimize.minimize cannot produce satisfactory results for my problem. Recently, I figured that CONOPT solver from GAMS is appropriate for my problem and following that I came across the suggestion for Pyomo last week as I was looking for some similar solver in Python. The problem I am trying to optimize is complex and not easy to setup in Pyomo for a new user, so I would appreciate if someone can help in setting up this problem using the model information and the description of its setup with scipy.optimize.minimize.

Problem:

The problem I am trying to solve involves one main equation, shown below in discretized form, that calls several other intermediate equations to update the model's parameters with time and return an output for each time. 

The number of parameters I want to optimize is 4 (a, b, and two other parameters that are present in the nested functions called from within the main equation), but each parameter is an array where the array size could vary depending on the complexity of the model. Below, I have described in brief how I setup this problem to solve with scipy.optimize.minimize.

Problem setup with scipy.optimize.minimize:

The objective function tries to minimize the squared difference between O_est and O_obs...where, O_est is the predicted array from fn_model() function that I am trying to match to the known array O_obs. I am interested in finding the optimized variable x_guess (a list of four variables) that minimizes the objective function as follow:

#=========== objective function that needs to be minimized
def min_fun_obj(x_guess):
    q_est = fn_main(x_guess) # q_est contains 83 elements
    fun_obj = np.sum([(q_o - q_e)**2 for q_o, q_e in zip(q_obs, q_est)])
    return fun_obj

The above minimization is dependent on function fn_model(), which calls 4 other nested functions, as shown:

#==============================================================================
# Main model function
#==============================================================================
def fn_model(x_guess):

    x1, x2, x3, x4 = x_guess[0]*np.ones((1)), x_guess[1]*np.ones((1, 1)), \
    x_guess[2]*np.ones((1)), x_guess[3]*np.ones((1))
    
    for i_p in xrange(n_p):
        for i_t in xrange(n_t-1):
            dt = t[i_t + 1] - t[i_t]
            #=============== call function fn1
            outp1 = fn1(x3[i_p])

            #=============== call function fn2
            O_est[i_t + 1, i_p] = fn2(x1[i_p], O_est[i_t, i_p], x2[:, i_p])
            #=============== call function fn3
            outp2 = fn3(data_est[i_t + 1, i_p])
            #=============== call function fn4
            outp3 = fn4(x3[i_p], dt, x4[i_p], x2[:, i_p], data_est[:, i_p])
           
    # outp1, outp2, outp3 are used within the model, but to keep it simple I am just showing them as outputs from the functions
    return O_est

I can minimize this problem (though unsatisfactorily) using scipy.optimize.minimize with bounds and constraints as shown below:

#=========== input for minimization
x_guess = [x[0], x[1], x[2], x[3]]) # initial guess
x_bounds = [(0, None), (0, 1), (0, 1), (0, None)]
cons = [{'type':'ineq', 'fun': lambda x_guess: x_guess[0]}, \
{'type':'ineq', 'fun': lambda x_guess: x_guess[1]}, \
{'type':'ineq', 'fun': lambda x_guess: (1 - np.sum(x_guess[3]))}]
#=========== call the scipy's minimization algorithm
x_optimized = optimize.minimize(min_fun_obj, x_guess, constraints=cons, bounds=x_bounds, method='SLSQP', jac = False, options={'disp': True})

[Siirola, John D]

This looks like a DAE parameter estimation problem.  My recommendation would be to NOT start with the discretized DAE and external functions, but rather to explicitly express the dynamic model in Pyomo using the pyomo.dae toolbox.  For examples, see the online documentation (https://software.sandia.gov/downloads/pub/pyomo/PyomoOnlineDocs.html#_dae_toolbox), and the parameter estimation example (https://github.com/Pyomo/pyomo/tree/master/examples/dae, look at Parameter_Estimation.py, run_Parameter_estimation.py, and data_set2.dat).

 

Depending on what your functions do, you may not even need to re-write them in Pyomo rules: if the function is a normal Python function that just performs simple algebra on the input variables (with no strange conditional logic or “weird” functions), passing in Pyomo variables may be sufficient to generate the necessary Pyomo expressions (that is, you can use your function as part of a constraint rule).

 

john

 

From: pyomo...@googlegroups.com [mailto:pyomo...@googlegroups.com] On Behalf Of Harpreet Singh
Sent: Wednesday, December 28, 2016 11:27 AM
To: Pyomo Forum <pyomo...@googlegroups.com>
Subject: [EXTERNAL] Help with setting up a complex optimization problem in pyomo

 

Hi. I have been looking for an optimization package for a month now after I figured scipy.optimize.minimize cannot produce satisfactory results for my problem. Recently, I figured that CONOPT solver from GAMS is appropriate for my problem and following that I came across the suggestion for Pyomo last week as I was looking for some similar solver in Python. The problem I am trying to optimize is complex and not easy to setup in Pyomo for a new user, so I would appreciate if someone can help in setting up this problem using the model information and the description of its setup with scipy.optimize.minimize.

 

Problem:

The problem I am trying to solve involves one main equation, shown below in discretized form, that calls several other intermediate equations to update the model's parameters with time and return an output for each time. 

 

--
You received this message because you are subscribed to the Google Groups "Pyomo Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Harpreet]

Thanks for the suggestion, John. I am doubtful that I would be able to use the original coupled ODEs because the parameter estimation problem from this model has been prescribed on the discretized form as presented in the original post. There are other intermediate functions that are either analytical expressions or a discretized equation and they perform the function of coupling the original ODEs...in my post earlier I am calling them as nested functions because they are intermediate functions that complete the overall model. 

Can I send you the code to help figure what problem I am trying to solve? I can send two codes for the problem setup to run using scipy and pyomo: the first code is setup to use scipy.optimize.minimize algorithm (it runs but the results are not satisfactory) and the second code is setup to use pyomo (this one is not running yet). Please let me know if I could send them by email.

Thanks.

[SH]

John, I took a look at the ODEs for the model and it seems that the discretized model has few parameters that come from some transformations and those parameters are not present in original ODEs. I would appreciate if any of you or someone else familiar with pyomo could help me figure out how to model my problem correctly.