# Python Example - Optimization

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

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.

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

$\displaystyle{ \sum_i \left(\frac{x_i-z_i}{z_i}\right)^2 = \sum_i \left(\frac{x_i}{z_i}-1\right)^2 }$

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!

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 ("Install pywin32:  https://pypi.org/project/pywin32/")
sys.exit(1)

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():
##
time.sleep(3)
del(Solver)
del(Tags)
del(Prj)