Python Example - Optimisation

From SysCAD Documentation
Jump to navigation Jump to search

Navigation: User Guide ➔ COM Automation ➔ Python Automation ➔ Example - Optimisation

Python Setup Python Examples
Installation &
Troubleshooting
Python
Utilities
Basic Usage & Scenarios Constrained FEM
(numpy|scipy|matplotlib)
Optimisation
(COM | scipy tools)
Model Testing
Framework
Dynamic
with GUI
Dynamic
External DLL
Programmatic
Model Generation
Importing data
to SysCAD

Related Links: Example Project - ChemApp with Python Optimisation


Overview

This example demonstrates using a number of tools to tune a steady state model (SysCAD project) based on data. This uses the Demo Precip Circuit in the distributed examples.

Full PSD Precipitation.png

We are provided with a product size distribution (for example from plant data) and we want to adjust some key parameters to tune the model to this data.

The two parameters we will adjust are the Agglomeration Rate and Cyclone Cut Size. Add the latter to the Trend Window so you can watch how these change as the optimisation proceeds.

In some ways this is an example of curve fitting - which is just another optimization problem. The curve fitting problem - find the "best fit" of a function with a small number of parameters to a set of data points, which means minimising the difference between the given data and the fit in some norm. Usually curve fitting involves evaluating a simple function to generate the fit. In this case, we generate the fit by running SysCAD while varying a small number of parameters.

Calculating the Objective Error

A key concept of optimisation is an objective function which is a target we want to minimise. In this case the target is to minimise the difference between the desired PSD and one our model calculates. We can calculate the objective function either in a PGM or in python - in an earlier version of this example we showed this in PGM, but we will now do the calculation in python. (A PGM listing is here, where the error is calculated and then fetched to the optimiser as a single value)

  • A more common target is minimisation of some "cost" associated with a particular configuration, or maximisation of "revenue". This would take into account total production (revenue) minus costs (for energy and raw materials).
  • In this example we calculate the objective function in python by fetching the PSD distribution using COM.
target = np.array([2.97869e-05, 3.93492e-05, 5.69551e-05, 8.58349e-05, 0.00013228, 0.00020766,
    0.00033343, 0.00055317, 0.00096398, 0.0018091, 0.0037611, 0.0087440, 0.022172,
    0.058889, 0.159535, 0.440231, 1.256861, 3.695465, 10.19094, 21.37749, 27.96972,
    21.47619, 9.873492, 2.837389, 0.542968, 0.073854, 0.0080610])

        psdist = np.array(sc.getTags(PSDTags))
        err = np.linalg.norm(psdist/target-1)

Python Optimisation

We will use the scipy optimisation module as well as numpy, so we import those. (These are part of the Anaconda distribution, otherwise pip install numpy/scipy)

import numpy as np 
import scipy.optimize as spo

The scipy.optimize module has a function fmin, which is the the Nelder-Mead downhill simplex algorithm

>>> help(spo.fmin)
Help on function fmin in module scipy.optimize.optimize:
fmin(func, x0, args=(), xtol=0.0001, ftol=0.0001, maxiter=None, maxfun=None, full_output=0, disp=1, retall=0, callback=None, initial_simplex=None)
   Minimize a function using the downhill simplex algorithm.   
   This algorithm only uses function values, not derivatives or second
   derivatives.

This takes a function (in python speak callable) (which is the objective function we want to minimise) and a starting point, plus a number of optional further parameters.

The callable func takes a single array argument func([x1, x2...]) , and evaluates the objective at these values, returning this result. In our case here, evaluation of the objective function means running SysCAD with the tags we want to adjust for opimisation set to these values, and return our calculated error. We have a list of optimising tag strings, and also need to fetch the outlet size distribution:

optTags =  ["PRECIP_CTRL.AggRateCorrection",  "CYC_001.d50Reqd (um)"]
PSDTags = ["P_011.Qo.Sz.I%d.FP.Al[OH]3(s) (%%)" % i for i in range(1, 28)]

Then our objective function is just

    def f(z):
        sc.setTags(optTags, z)
        sc.run()
        psdist = np.array(sc.getTags(PSDTags))
        err = np.linalg.norm(psdist/target-1)
        sc["PlantModel.PUV1.Value"]=err        ## We store this back in SysCAD so we can plot the error in the trend window
        return err

Having defined this, then we just make a call to fmin with some suitable starting values (in the code, just fetch the existing tags)

res = spo.fmin(f, ystart)

Thats all there is to it!

Pyopt01.png

The desired PSD distribution was actually determined by running the model with agglomeration correction 0.4 and cut size 70 microns, and these valued are recovered by the optimization.

Be sure to call doClose() after finishing to clean up COM objects, otherwise you end up with zombie SysCAD instances. (If strange things are happening use Process Explorer or Task Manager to remove these)

Optimization terminated successfully.
         Current function value: 0.000053
         Iterations: 39
         Function evaluations: 74
[ 0.40000361 69.99949443]
Full Python Code      
#! python3 -i
## Optimisation with COM and python.

import syscadif
import numpy as np
import scipy.optimize as spo

ScdPrj = r'C:\SysCAD139\Examples\10 Alumina\PSD Precipitation Circuit Example.spf\Project.spj' 
optTags =  ["PRECIP_CTRL.AggRateCorrection",  "CYC_001.d50Reqd (um)"]
PSDTags = ["P_011.Qo.Sz.I%d.FP.Al[OH]3(s) (%%)" % i for i in range(1, 28)]

target = np.array([2.97869e-05, 3.93492e-05, 5.69551e-05, 8.58349e-05, 0.00013228, 0.00020766,
0.00033343, 0.00055317, 0.00096398, 0.0018091, 0.0037611, 0.0087440, 0.022172,
0.058889, 0.159535, 0.440231, 1.256861, 3.695465, 10.19094, 21.37749, 27.96972,
21.47619, 9.873492, 2.837389, 0.542968, 0.073854, 0.0080610])

sc = syscadif.SysCADCom()
sc.OpenProject(ScdPrj)
 
def main():
    ## start from saved solution
    ystart = np.array(sc.getTags(optTags))
    print (ystart)
 
    def f(z):
        sc.setTags(optTags, z)
        sc.run()
        psdist = np.array(sc.getTags(PSDTags))
        err = np.linalg.norm(psdist/target-1)
        sc["PlantModel.PUV1.Value"]=err   ## We store this back in SysCAD so we can plot the error in the trend window
        return err
 
    res = spo.fmin(f, ystart, xtol=.001, ftol=.05)
    print (res)

main()
input('Enter to finish')
sc.Close()

Remarks

  • You could also do optimisation based on tuning external parameters instead of SysCAD tags, for example in PHREEQC.
  • The objective function could also be calculated in PGM, e.g.:
Objective Function in PGM      
;--- variable declarations ---
PageLabel(Precip)
    TextLabel(,"Recycle mass flows")
    REAL 	FS_SolidRecovery*("Frac", "%")
    REAL 	FS_Solids_Moisture*("Frac", "%")
    REAL 	CS_SolidRecovery*("Frac", "%")
    REAL 	CS_Solids_Moisture*("Frac", "%")
    TextLabel(,"Global Agglomeration correction")
    Long	i, NumberofTanks*<<14>>
    REAL 	AggRateCorrection*
    TextLabel(,"Yield Calculations")
    REAL 	Production@("Qm", "t/d")
    REAL 	Yield@("Conc", "g/L")
    REAL 	ProductD50@("L", "um")
    TextLabel(,"Check SysCAD Version", "Function IntStr0 available in SysCAD9.3 Build136.19015")
    Checkbox BuildOK@
 
;; Constant Array with desired product size distribution	
Array z = {2.97869e-05, 3.93492e-05, 5.69551e-05, 8.58349e-05, 0.00013228, 0.00020766,
0.00033343, 0.00055317, 0.00096398, 0.0018091, 0.0037611, 0.0087440, 0.022172,
0.058889, 0.159535, 0.440231, 1.256861, 3.695465, 10.19094, 21.37749, 27.96972,
21.47619, 9.873492, 2.837389, 0.542968, 0.073854, 0.0080610 }
    
;; Work array 
Array x
PageLabel(Analysis)
TextLabel(, "Error difference between product PSD and desired PSD") 
real PSDError@
    	
Sub InitialiseSolution()
    ["X_003.Separ.SolidsToUFReqd (%)"]    = FS_SolidRecovery
    ["X_003.Separ.UFSolidFracReqd (%)"] = 100 - FS_Solids_Moisture
    ["X_001.Separ.SolidsToUFReqd (%)"]    = CS_SolidRecovery
    ["X_001.Separ.UFSolidFracReqd (%)"] = 100 - CS_Solids_Moisture
    
    i = 1
    While i <= NumberofTanks
    	[Concatenate("PC_", IntStr0(i, 3), ".Agglom.Rate.Correction")] = AggRateCorrection
    	i = i + 1
    EndWhile
    x.SetLen(27)
EndSub

    Production =  ["P_011.Qo.SQm (t/d)"] 
    Yield = Production/["P_012.Qv (kL/d)"]*1000.*102./156.
    ProductD50 = ["P_011.Qo.Sz.d50.Al[OH]3(s) (um)"]
    i = 1
    while (i<=27) 
       x[i-1] = [Concatenate("P_011.Qo.Sz.I", IntToStr(i), ".FP.Al[OH]3(s) (%)")]
       i = i+1
    endwhile
    x.Div(z)
    x.Offset(-1)
    PSDError = x.inner(x)
        
$ ; --- end of file ---