Python Example - Optimization

From SysCAD Documentation
Jump to navigation Jump to search

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

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: Full PGM code for Python Optimization Example


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.

Start by making a copy of this project (Save Version) since we will be modifying the PGM, and work with the copy.

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 optimization 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 minimizing 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 optimization is an objective function which is a target we want to minimize. In this case the target is to minimize 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 this case we will do it in the PGM, using some of the built-in array manipulation class functionality.

  • A more common target is minimization of some "cost" associated with a particular configuration, or maximization of "revenue". This would take into account total production (revenue) minus costs (for energy and raw materials).

Edit the PGM for the model and add the following lines to the declaration section:

;; 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 [email protected]

In the InitialiseSolution() subroutine, set the x-array length to 27

	x.SetLen(27)

Now in the logic executed at every step, get the PSD distribution and calculate the difference:

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)

Here we are calculating

[math]\displaystyle{ \sum_i \left(\frac{x_i-z_i}{z_i}\right)^2 = \sum_i \left(\frac{x_i}{z_i}-1\right)^2 }[/math]

using array operations - divide x by z, subtract 1, then calculate the inner product with itsself (which is just the sum of squares)

Full PGM listing is here

  • In this case we are calculating the objective function in SysCAD, we could also calculate the objective function in python by fetching the PSD distribution using COM.

Python Optimization

We will use the scipy optimization 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 minimize) 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 opimization set to these values, and return our calculated error. We have a list of optimizing tag strings, and an objective tag string:

optTags =  ["PRECIP_CTRL.AggRateCorrection",  "CYC_001.d50Reqd (um)"]
objTag = "PRECIP_CTRL.PSDError"

Then our objective function is just

def f(z):
  for  v, theTag in zip(z, optTags): ## run through pairs of tag-value
    Tags.SetTagValue(theTag, v)     
  Run()                              ##   
  return Tags.TagValue(objTag)       ##  fetch the objective tag 

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.000000
        Iterations: 39
        Function evaluations: 74
[ 0.40000361 69.99949443]
>>> doClose()

Code

## Optimization with COM and python.
import sys
import time
import numpy as np
import scipy.optimize as spo
try:
  import win32com.client as wc    ## PyWin32 COM Client
except:
  print ("Win32 python extensions not found")
  print ("Install pywin32:  https://pypi.org/project/pywin32/")
  sys.exit(1)

SysCAD = wc.DispatchEx("SysCADSimulator93.Application")  ## Fire up SysCAD
ScdDir = r'C:\SysCAD138'
ScdPrj = r'\Examples\10 Alumina\PSD Precipitation Circuit Example(Optimizer Test)-01.spf\Project.spj'
 
Prj = SysCAD.OpenProject(ScdDir+ScdPrj)  ## Open project 
Tags = Prj.Tags
Solver = Prj.Solver
ProBal = Solver.ProBal
 
optTags =  ["PRECIP_CTRL.AggRateCorrection",  "CYC_001.d50Reqd (um)"]
objTag = "PRECIP_CTRL.PSDError"
 
def Run():
    ##Start Probal and wait until it is solved
    global iter
    print(iter, end="")
    ProBal.Start()                        ## Start ProBal
    while True:                           ## Wait until solved 
       print (".", end="")                 
       time.sleep(.1)
       if ProBal.IsStopped:
           iter+=1
           print()
           break
 
def doGo():
    ## start from saved solution
    ystart = np.zeros(len(optTags))
    for i, tag in enumerate(optTags): 
        ystart[i] = Tags.TagValue(tag)
    print (ystart)
 
    def f(z):
        for v, theTag in zip(z, optTags):
            Tags.SetTagValue(theTag, v)
        Run()
        return Tags.TagValue(objTag)
 
    res = spo.fmin(f, ystart, xtol=.001, ftol=.05)
    print (res)
 
def main():
    global iter
    iter = 0
    doGo()
    
def doClose():
    global SysCAD, Solver, Tags, Prj
    SysCAD.CloseProject(False)
    ##
    time.sleep(3)
    del(Solver)
    del(Tags)
    del(Prj)
    del(SysCAD)
 
   
main()

Remarks

  • You could also do optimization based on tuning external parameters instead of SysCAD tags, for example in PHREEQC