Embedded Python in Modelica® am Beispiel einer Optimierung mit Neben- und Randbedingung in SimulationX

Die folgenden Notizen zeigen wie eine in Python implementierte Funktion bzw. ein Python Modul aus Modelica® aufgerufen werden kann. Dies erfolgt über Pure Embedding (siehe https://docs.python.org/2/extending/embedding.html#pure-embedding) und einem C++-Wrapper, der die Python Funktion als eine externe Funktion kapselt. Implementierung und Test erfolgten mit SimulationX.

Example for optimization with scipy.optimize using SLSQP method

The original problem formulation is an example by Matlab in Chemical Engineering at CMU using fmincon in MATLAB

http://matlab.cheme.cmu.edu/2011/12/25/finding-equilibrium-composition-by-direct-minimization-of-gibbs-free-energy-on-mole-numbers/

Library imports

In [1]:
from numpy import *
from scipy.optimize import minimize

Definition of constants for later test

In [2]:
R=8.3144598
Aeq=array([[1,0,1,1,0,0,0,2,3,4],[4,2,0,0,2,0,0,6,8,10],[0,1,1,2,0,0,2,0,0,0],[0,0,0,0,0,2,0,0,0,0]])
B=array([[0.230769230769231],[0.615384615384615],[0.115384615384615],[0.0384615384615385]])
beq=B[:,0]
T=450+273.15
p=1.01325e5
p0=101325
z0=ones(10)*0.1
G0=array([-218418.269410093,-385823.782728340,-259833.834235695,-557093.447589739,-100798.146587938,
          -144836.255076941,-154900.925575189,-264145.211175915,-320966.890140181,-377373.644866805])

Boundary for $z$

In [3]:
bounds=vstack((ones(10)*1e-15,ones(10)*1)).transpose()
bounds
Out[3]:
array([[  1.00000000e-15,   1.00000000e+00],
       [  1.00000000e-15,   1.00000000e+00],
       [  1.00000000e-15,   1.00000000e+00],
       [  1.00000000e-15,   1.00000000e+00],
       [  1.00000000e-15,   1.00000000e+00],
       [  1.00000000e-15,   1.00000000e+00],
       [  1.00000000e-15,   1.00000000e+00],
       [  1.00000000e-15,   1.00000000e+00],
       [  1.00000000e-15,   1.00000000e+00],
       [  1.00000000e-15,   1.00000000e+00]])

Definition of objective function $G$ and its Jacobian $\frac{\partial G}{\partial \underline {z}}$

The objective function is: $$\begin{align} G_{\min}&=\sum\limits_i z_i \left(\frac{G_{0,i}}{R T} + \log\left(\frac{z_i}{\sum_j{z_j}}\frac{p}{p_0}\right)\right)\\ &=\sum\limits_i \left( z_i\frac{G_{0,i}}{R T} + z_i \log\left(z_i\frac{p}{p_0}\right)-z_i\log\left(\sum_j z_j\right)\right) \end{align}$$

The gradient in the direction of $k$ ($\frac{\partial}{\partial z_k}$): $$\begin{align} \frac{\partial}{\partial z_k} G_{\min}&= \frac{\partial}{\partial z_k} \sum\limits_i \left( z_i\frac{G_{0,i}}{R T} + z_i \log\left(z_i\frac{p}{p_0}\right)-z_i\log\left(\sum_j z_j\right)\right)\\ &= -\sum\limits_{i, i\neq k}\left(\frac{z_i}{\sum_j z_j}\right) + \frac{G_{0,k}}{R\,T}+\log\left(z_k\frac{p}{p_0}\right)+1-\frac{z_k}{\sum_j z_j}-\log\sum_j z_j\\ &=\frac{G_{0,k}}{R\,T}+\log\left(\frac{z_k}{\sum_j z_j}\frac{p}{p_0}\right) \end{align}$$ (with $\sum\limits_{i, i\neq k}\left(\frac{z_i}{\sum_j z_j}\right)+\frac{z_k}{\sum_j z_j}=\frac{\sum_i z_i}{\sum_j z_j}=1$)

In [4]:
def Gmin(z, G0, R, T, p, p0):
    z[z<=0]=1e-20
    return sum(z*(G0/R/T + log(z*p/p0)-log(sum(z))))

def JacGmin(z, G0, R, T, p, p0):
    dz=0*z
    z[z<=0]=1e-20
    for k,tz in enumerate(z):
        dz[k]=G0[k]/(R*T) + log(z[k]*p/p0)-log(sum(z))
    return dz[:]
    

Definition of constraint $\mathbf A \underline z - \underline b=0$ and its Jacobian which is $A$

In [5]:
def AzB(z, A, b):
    return dot(A,z)-b

def JacAzB(z, A, b):
    return A

Callback function (not used)

In [6]:
def print_z(z):
    print z

Call of minimize

In [7]:
def call_optimization(z0,G0,R,T,p,p0,Aeq,beq):
    cons1 = {'type': 'eq', 
         'fun': AzB, 
         'jac': JacAzB, 
         'args': (Aeq, beq)}
    res = minimize(Gmin,
                z0,  
                args=(G0, R, T, p, p0),
                method='SLSQP',
                jac=JacGmin,
                hess=None,
                hessp=None,
                bounds=bounds,
                constraints=cons1,
                tol=1e-10,
                callback=None, #print_z,
                options={'disp': False, 'iprint': 1, 'maxiter': 100, 'ftol': 1e-10})
    #print res['message']
    #print res
    z=res['x']
    #print 'constraint:'
    #print AzB(z, Aeq, beq)
    return z[:]

Only for test of call_optimization

In [8]:
z=call_optimization(z0,G0,R,T,p,p0,Aeq,beq)
Gmin(z, G0, R, T, p, p0)
Out[8]:
-11.551364414178844
In [10]:
z = array([0.151227518198636,0.000488474630907019,0.0440616190347329,0.0354172608594877,0.00465456614012380,
           0.0192307692307692,1.82984720770906e-21,3.13608665175083e-05,3.69811129496762e-08,
           1.32526281497577e-21])
Gmin(z, G0, R, T, p, p0)
Out[10]:
-11.551364413165267

Optional install python module using setup.py

# -*- coding: utf-8-*-
from distutils.core import setup

setup(name='optimization_Gmin',
    version='1.0',
    description = 'optimization of a certain objective',
    author='KAO, ESI ITI GmbH, Dresden Germany',
    url='http://www.simulationx.com',
    py_modules=['optimization_Gmin'])

Notes

  • set includes for <Python.h> and <numpy\arrayobject.h>
  • link against npymath.lib and python27.lib
  • import_array() must be called within the .dll only once. Therefore a static bool variable might be used. This shouldn't be called within DLL_PROCESS_ATTACH but in a separately defined method.
  • Py_Initialize() must be called iff !Py_IsInitialized()
  • the pointers PyObject* which point on Python modules of functions should be static. When the .dll is detached (DLL_PROCESS_DETATCH) the pointers should be dereferences and set to NULL. Loading of modules can be done in the exported function but only if the static pointers (PyObject*) are NULL.
  • The Python Modul is located in the SimulationX working directory and is therein translated into .pyc.
  • During a simulation DLL_PROCESS_ATTACH and DETATCH can be called multiple times. Without additional statements the PyObject* which point at modules or funktions keep their value and would be valid. After a reset they are NULL, because the DLL is reloaded.
  • PyObject* which are passed as arguments to functions using a tuple must not separately be dereferenced using Py_DECREF(), this is done already by dereferencing the tuple.
  • Achanges in pyconfig.h to be able to builöd debug version of the .dll:
    • Uncomment #define Py_DEBUG i.e.: //#define Py_DEBUG
    • remove python27_d.lib or change into python27.lib
# ifdef _DEBUG 
#     pragma comment(lib,"python27_d.lib") 
# else 
#     pragma comment(lib,"python27.lib") 
# endif /* _DEBUG */

in

# ifdef _DEBUG 
#     pragma comment(lib,"python27.lib") 
# else 
#     pragma comment(lib,"python27.lib") 
# endif /* _DEBUG */

Implementation

Existing Package in the Buildings package

Link: http://simulationresearch.lbl.gov/modelica/

Buildings.Utilities.IO.Python27.Real_Real

Information Sources

Python Online Documentation: Chapter 5. Embedding Python in Another Application https://docs.python.org/2/extending/embedding.html

External Function definition in Modelica

function MinimizeG "Example for external function z0,G0,R,T,p,p0,A,b"
    input Real z0[:];
    input Real G0[:];
    input Real R;
    input Real T;
    input Real p;
    input Real p0;
    input Real A[:,:] "Input Parameter";
    input Real b[:] "Input Parameter";
    output Real z[size(z0,1)] "Output Parameter";
    external "C" dllPyWrapperMinimize(z0,size(z0,1),G0,size(G0,1),R,T,p,p0,A,size(A,1),size(A,2),b,size(b,1),z,size(z,1))
        annotation(__iti_dll="PyWrapperMinimize.dll");
end MinimizeG;

Function call within Modelica

G_min = sum(z.*(G0/R/T + log(z/sum(z)*p/p0)));
z = ExternalCodeTraining.MinimizeG(z0,G0,R,T,p,p0,Aeq,vector(Beq));
psi_eq=z/sum(z);

C++ Wrapper Source

// PyWrapperMinimize.cpp : Defines entry point into dll.
//

#include "stdafx.h"
#include <Python.h>
#include <numpy\arrayobject.h>
//#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION

#ifdef _MANAGED
#pragma managed(push, off)
#endif

static bool bInitialized;

typedef void TRACE_FUNCTION(int, int, const char*, const char*);
static TRACE_FUNCTION* trace = NULL;
static HINSTANCE hInstGlobalHlp = NULL;

static int GlobalCounter = 0;

static PyObject *pName = NULL;
static PyObject *pModule = NULL;
static PyObject *pDict = NULL;
static PyObject *pFunc = NULL;

void traceFunction(int kind, const char* title, char* msg);

void MyPyInit()
{
    if (!bInitialized) {
        import_array();
        bInitialized = true;
    }
    return;
}

BOOL APIENTRY DllMain(HMODULE hModule,
    DWORD  ul_reason_for_call,
    LPVOID lpReserved
    )
{
    switch (ul_reason_for_call)
    {
    case DLL_PROCESS_ATTACH:
    {
            hInstGlobalHlp = LoadLibrary(L"globalhlp.dll");
            if (hInstGlobalHlp)
            {
                trace = (TRACE_FUNCTION*)GetProcAddress(hInstGlobalHlp, "OutputTraceMessageA");
            }
            break;
    }
    case DLL_THREAD_ATTACH:
    case DLL_THREAD_DETACH:
    case DLL_PROCESS_DETACH:
    {

        Py_XDECREF(pFunc);
        Py_XDECREF(pModule);
        pFunc = NULL;
        pModule = NULL;
        if (hInstGlobalHlp)
        {
            FreeLibrary(hInstGlobalHlp);
            hInstGlobalHlp = NULL;
            trace = NULL;
        }
        break;
    }
    }
    return TRUE;
}


/* This is the function that is called from SimulationX */
extern "C" __declspec(dllexport) void dllPyWrapperMinimize( double* z0, size_t rows_z0,
                            double* G0, size_t rows_G0,
                            double R,
                            double T,
                            double p,
                            double p0,
                            double* A, size_t rows_A, size_t cols_A,  
                            double* b, size_t rows_b,
                            double* z, size_t rows_z
                           )
{

        PyObject *pArgs;
        PyObject* retVal = NULL;

        if (!Py_IsInitialized())
            Py_Initialize();

        MyPyInit();
        if (!pModule && !pName) {
            pName = PyString_FromString("optimization_Gmin");
        }
        if (!pModule) {
            // Here the Module (.py) is loaded
            pModule = PyImport_Import(pName);
        }
        Py_XDECREF(pName);
        pName = NULL;


        int succsess=1;
        // Error checking of pName left out
        // fprintf(stderr, "Module to load %s\n", argv[2]);

        /*
        PyObject* PyA = PyArray_NewFromDescr(&PyArray_Type,
                            PyArray_DescrFromType(NPY_DOUBLE),
                            2, // matrix
                            py_shapeA, NULL,
                            (void *)0,
                            NPY_ARRAY_OWNDATA | NPY_ARRAY_CARRAY,
                            NULL);
        */
        /*
        double* p_dbl_A = (double*)PyArray_DATA(PyA);
        for (int i = 0; i < py_shapeA[1]; i++) {
            for (int k = 0; k < py_shapeA[0]; k++) {
            p_dbl_A[k + i*py_shapeA[0]] = A[i + k*py_shapeA[1]]; // careful
            }
        }
        */

        if (pModule != NULL) {
            if (pFunc == NULL){
                // Here the function is loaded from the module
                pFunc = PyObject_GetAttrString(pModule, "call_optimization");
            }
            /* pFunc is a new reference */

            if (pFunc && PyCallable_Check(pFunc)) {
                pArgs = PyTuple_New(8);

                npy_intp py_shape_z0 = rows_z0;  // shape of matrix
                PyObject* Py_z0 = PyArray_SimpleNewFromData(1, &py_shape_z0, NPY_DOUBLE, z0);

                npy_intp py_shape_G0 = rows_G0;  // shape of matrix
                PyObject* Py_G0 = PyArray_SimpleNewFromData(1, &py_shape_G0, NPY_DOUBLE, G0);

                PyObject* Py_R = PyFloat_FromDouble(R);
                PyObject* Py_T = PyFloat_FromDouble(T);
                PyObject* Py_p = PyFloat_FromDouble(p);
                PyObject* Py_p0 = PyFloat_FromDouble(p0);

                npy_intp py_shape_A[2];  // shape of matrix
                py_shape_A[0] = rows_A;
                py_shape_A[1] = cols_A;
                PyObject* Py_A = PyArray_SimpleNewFromData(2, py_shape_A, NPY_DOUBLE, A);

                npy_intp py_shape_b = rows_b;  // shape of matrix
                PyObject* Py_b = PyArray_SimpleNewFromData(1, &py_shape_b, NPY_DOUBLE, b);

                /*for (i = 0; i < argc - 3; ++i) {
                pValue = PyInt_FromLong(atoi(argv[i + 3]));
                if (!pValue) {
                Py_DECREF(pArgs);
                Py_DECREF(pModule);
                fprintf(stderr, "Cannot convert argument\n");
                return 1;
                }
                // pValue reference stolen here:
                PyTuple_SetItem(pArgs, i, pValue);
                }*/
                succsess = PyTuple_SetItem(pArgs, 0, (PyObject*)Py_z0);
                succsess = PyTuple_SetItem(pArgs, 1, (PyObject*)Py_G0);
                succsess = PyTuple_SetItem(pArgs, 2, (PyObject*)Py_R);
                succsess = PyTuple_SetItem(pArgs, 3, (PyObject*)Py_T);
                succsess = PyTuple_SetItem(pArgs, 4, (PyObject*)Py_p);
                succsess = PyTuple_SetItem(pArgs, 5, (PyObject*)Py_p0);
                succsess = PyTuple_SetItem(pArgs, 6, (PyObject*)Py_A);
                succsess = PyTuple_SetItem(pArgs, 7, (PyObject*)Py_b);

                ////////////////////////////////////////////////
                // That's the function call to Python 
                ////////////////////////////////////////////////
                retVal = PyObject_CallObject(pFunc, pArgs);

                Py_DECREF(pArgs);

                /* This must not be done, its done Py_DECREF(pArgs). If its done before 2nd evaluation will crash 
                Py_XDECREF(Py_z0);
                Py_XDECREF(Py_G0);
                Py_XDECREF(Py_R);
                Py_XDECREF(Py_T);
                Py_XDECREF(Py_p);
                Py_XDECREF(Py_p0);
                Py_XDECREF(Py_A);
                Py_XDECREF(Py_b);
                */


                npy_intp py_shape_z = rows_z;
                if (retVal != NULL) {
                    double* po = (double*)PyArray_DATA(retVal);
                    for (int i = 0; i < rows_z; i++) {
                        z[i] = po[i];
                    }
                    Py_DECREF(retVal);
                }
                else {
                    PyErr_Print();
                    fprintf(stderr, "Call failed\n");
                    return;
                }

            }
            else {
                if (PyErr_Occurred())
                    PyErr_Print();
            }

        }
        else {
            PyErr_Print();
            return;
        }

        return;

}

#ifdef _MANAGED
#pragma managed(pop)
#endif

Outlook

Same function using embedding

  • Scilab
  • Julia language
In [ ]: