Python Example - Optimisation

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

Python Setup	Python Examples	Python Script			Optimisation		Visulisation	Python GUI	Tags and Data
Installation & Troubleshooting	List of Examples	Simple Script (pywin32)	SysCAD COM Python Class	Automated Model Testing	Constrained FEM (numpy\|matplotlib)	Optimisation (numpy\|scipy)	Adding Plots (numpy\|matplotlib)	Dynamic with GUI	Accessing Data (sqlite3\|pandas)

Overview

This example shows how to use a range of tools to tune a steady-state SysCAD model using real-world data. It’s based on the PSD Precip with Oxalate co-precip Example project.

In many ways, this is a form of curve fitting, which is essentially an optimisation problem. The goal is to find the best fit of a model (with a small number of adjustable parameters) to a set of data points by minimising the difference between the model output and the actual data. Unlike traditional curve fitting, which often uses simple mathematical functions, this approach involves running SysCAD simulations while varying parameters.

We’re working with a product size distribution, typically derived from plant data, and our goal is to adjust key model parameters to better align the simulation output with this observed data.

The two parameters we’ll be adjusting are:

Agglomeration Rate
Cyclone Cut Size

In this example, the project was previously solved using a specific set of parameters. The optimisation process is used here to recover those same values. As a result, the optimised output matches the target data exactly as shown in the plot below.

Calculating the Objective Error

A key concept in optimisation is the objective function — a target we aim to minimise. In this case, the objective is to reduce the difference between the target PSD and the PSD predicted by the model.

Earlier versions of this example calculated the objective function using a PGM (see Objective Function in PGM). However, in this version, the calculation is performed in Python by retrieving the PSD data via the COM interface.

In broader applications, the objective function might represent a cost to minimise or a revenue to maximise—taking into account factors such as production output, energy consumption, and raw material usage. For this example, though, the focus remains on minimising the PSD error.

target = np.array([1.87e-05,3.86e-05,7.88e-05,0.0001574,0.00030658,0.00059523,
                   0.0012118,0.0026943,0.0064393,0.015705,0.03767,0.08799,0.20063,
                   0.44764,0.97731,2.0864,4.349,8.7497,16.047,23.145,22.01,13.522,
                   5.822,1.8915,0.48233,0.098148,0.0184])

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

Python Optimisation

## Optimisation with COM and python.
import syscadif
import numpy as np
import scipy.optimize as spo

# Create a SysCAD COM object to interact with the simulation
sc = syscadif.SysCADCom()
ScdDir = r'C:\SysCAD%d' % sc.build   ## Installation folder
ScdPrj = r'\Examples\10 Alumina\PSD Precip with Oxalate co-precip Example.spf\Project.spj' 

optTags =  ["PRECIP_CTRL.AggRateCorrection",  "CYC_001.d50Reqd (um)"] ## Tag to be optimised:
PSDTags = ["P_044.Qo.Sz.I%d.FP.Al[OH]3(s) (%%)" % i for i in range(1, 28)]

## Objective Function
target = np.array([1.87e-05,3.86e-05,7.88e-05,0.0001574,0.00030658,0.00059523,
                   0.0012118,0.0026943,0.0064393,0.015705,0.03767,0.08799,0.20063,
                   0.44764,0.97731,2.0864,4.349,8.7497,16.047,23.145,22.01,13.522,
                   5.822,1.8915,0.48233,0.098148,0.0184])

sc.OpenProject(ScdDir+ScdPrj)  # Open the SysCAD project
 
def main():
    ## Get initial values for the optimisation parameters from the simulation
    ystart = np.array(sc.getTags(optTags))
    print (ystart)
    ##Objective function to minimise: calculates error between simulated and target PSD
    def f(z):
        sc.setTags(optTags, z)                  # Update simulation parameters
        sc.run()                                # Run the simulation
        psdist = np.array(sc.getTags(PSDTags))  # Get PSD results
        err = np.linalg.norm(psdist/target-1)   # Calculate relative error
        sc["PlantModel.PUV1.Value"]=err         # Write error back to SysCAD for monitoring
        return err
    ## Run optimisation using Nelder-Mead method
    res = spo.fmin(f, ystart, xtol=.001, ftol=.05)
    print (res)

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

The scipy.optimize module includes a function called fmin, which uses the Nelder-Mead downhill simplex algorithm. It’s designed to minimise an objective function and requires:

A callable (a Python function) that represents the objective we want to minimise.
A starting point for the optimisation.
A few optional parameters.

The function you pass in should take a single array argument, like func([x1, x2, ...]), and return the result of evaluating the objective. In our case, that means running SysCAD with the optimisation tags set to those values, then returning the calculated error.

We’ve got a list of tag strings that we’re optimising, and we also need to retrieve the outlet size distribution.

Once that’s set up, we just call fmin with some sensible starting values—usually by grabbing the current tag values from the model. That’s all there is to it!

The target particle size distribution (PSD) was originally determined by running the model with an agglomeration correction of 3 and a cut size of 70 microns. These values are successfully recovered through the optimisation process.

[  4. 110.]  
Optimisation terminated successfully.
         Current function value: 0.002574
         Iterations: 44
         Function evaluations: 84
[ 2.99669813 70.00395296]

Adding the variables to the Trend Window in SysCAD will provide visualisation on how values change during the optimisation process.

Remarks

You could also do optimisation based on tuning external parameters instead of SysCAD tags, for example in PHREEQC. (See Discussion: Ternary Saturation Diagrams using PHREEQC)
The objective function could also be calculated in PGM, e.g.:

Objective Function in PGM

;--- variable declarations ---
PageLabel "Precip"
    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")

;; Constant Array with desired product size distribution	
Array z = {1.87e-05,3.86e-05,7.88e-05,0.0001574,0.00030658,0.00059523,
           0.0012118,0.0026943,0.0064393,0.015705,0.03767,0.08799,0.20063,
           0.44764,0.97731,2.0864,4.349,8.7497,16.047,23.145,22.01,13.522,
           5.822,1.8915,0.48233,0.098148,0.0184 }
    
;; Work array 
Array x
TextLabel , "Error difference between product PSD and desired PSD" 
real PSDError@
    	
Sub InitialiseSolution()
    i = 1
    While i <= NumberofTanks
    	[Concatenate("PC_", IntStr0(i, 3), ".Agglom.Rate.Correction")] = AggRateCorrection
    	i = i + 1
    EndWhile
    x.SetLen(27)
EndSub

    Production =  ["P_044.Qo.SQm (t/d)"] 
    Yield = Production/["P_100.Qv (kL/d)"]*1000.*102./156.
    ProductD50 = ["P_044.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 ---

Python Example - Optimisation

Contents

Overview

Calculating the Objective Error

Python Optimisation

Remarks

Navigation menu

Python Example - Optimisation

Overview

Calculating the Objective Error

Python Optimisation

Remarks

Navigation menu

Search