Python Example - Adding Visualisation with Plots

From SysCAD Documentation
(Redirected from Python Utilities)
Jump to navigation Jump to search

Navigation: User Guide ➔ COM Automation ➔ Python Automation ➔ Python Example - Adding Visualisation with Plots

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)

Python has lots of powerful libraries for data manipulation, plotting, and analysis. These are some examples of what can be done

Adding visualisation with Plots

As we run scenarios, it results are best presented using plots. Here we will discuss how to add visualisation to some example projects already set up to run scenarios.

Gold Example Project

The SysCADCom class has a useful function to run a series of scenarios and collect the values of any number of tags.
The project was run across a range of feed rates. From the resulting plot, it is evident that the current setup is not suitable for higher flow rates. To accommodate these higher flow rates, a larger tank volume is required. You may consider re-running the scenario with an increased tank volume in the CSTR configuration and reviewing the updated results.

import syscadif
import matplotlib.pyplot as plt

sc = syscadif.SysCADCom()
ScdDir = r'C:\SysCAD%d' % sc.build
ScdPrj = r'\Examples\25 Gold\Demo Gold Project.spf\Project.spj'
 
sc.OpenProject(ScdDir+ScdPrj)  ## Open project 

X = list(range(40, 300, 50))  ## for Python3 since since range() is now a generator and we need to reuse X
Y = sc.RunScenarios(X, "SLURRY_IN.QmReqd (t/h)", ["P_001.Qo.QEl.Au (oz/d)", "P_016.Qo.QEl.Au (oz/d)", "P_007.Qo.QEl.Au (oz/d)"])
Z = [y[0] for y in Y]  # flatten
Z1 = [y[1] for y in Y]
Z2 = [y[2] for y in Y]
plt.plot(X, Z, label="Feed")
plt.plot(X, Z, "ro")
plt.plot(X, Z1, label = "LoadedCarbon")
plt.plot(X, Z2, label = "Tails")
plt.xlabel("Slurry Feed Rate (t/h)")
plt.ylabel("Gold Flowrate (oz/d)")
plt.legend(loc=2)
plt.show()
 
sc.CloseProject()
sc.Close()
Comexample1.png

Washer Example

Here is the Washer Example Project showing the NaOH concentration profiles.

Python Code      
import numpy as np
import matplotlib.pyplot as plt

import sys
sys.path.append(r"C:\SysCAD139\BaseFilesUser\Scripts")   ## syscadif.py is saved in a subdirectory of BaseFilesUser
import syscadif

sc = syscadif.SysCADCom()
ScdDir = r'C:\SysCAD%d' % sc.build
ScdPrj = r'\Examples\03 UnitModels\Counter Current Washer Example.spf\Project.spj'

WashTags  = ["CCW_CONTROL.CCW{}_Online",
             "WASHER_CONTROL.CCD00{}_Online",
             "WASHER_CONTROL.CCD10{}_Online",
             "TH_CONTROL.TH{}_Online"]
OutTags = ["CCW_CONTROL.UF{}_NaOHConc (g/L)",
           "WASHER_CONTROL.UF00{}_NaOHConc (g/L)",
           "WASHER_CONTROL.UF10{}_NaOHConc (g/L)", 
           "TH_CONTROL.UF{}_NaOHConc (g/L)"]

def setStages(ion):
    '''Set all the on/off states for washers'''
    for i in range(2, 7):
        for tg in WashTags:
            sc[tg.format(i)] = (i<ion)

def getData():
    return np.array([sc[tg.format(i)] for tg in OutTags for i in range(1, 7)])
    
def t0():
    plt.ion()
    fig, axs = plt.subplots(1, 4, sharey=True, tight_layout=True)

    for i in range(2, 7):
        setStages(i)
        sc.run()
        dd = getData().reshape((4, 6))
        for ax, dat in zip(axs, dd):
            ax.plot(dat, label=f"{i} Stages")
            ax.plot(dat, "ko")   
    for ax in axs: ax.legend()
    for ax, ss in zip(axs, ["CC Washer", "Washer (Bypass to OF)", "Washer (Bypass to UF)", "Thickener"]):
        ax.set_title(ss)
    axs[0].set_ylabel("NaOH Conc (gpl)")
    plt.suptitle("NaOH Profiles")
    plt.savefig("fig4.png")
    plt.show()

sc.OpenProject(ScdDir + ScdPrj)
sc.run()

t0()

sc.CloseProject()
sc.Close()

WasherExample.png

Sensitivity Analysis Example

Here is the Sensitivity Analysis example done with python.

We could run this as a set of scenarios, but for fun we use a numpy meshgrid to create the points at which we want to evaluate the model. We also create a vectorized function which is used to evaluate the model at each of the points. We use the resulting 2D array to plot a 3D graphic. (See left images below.)

Python Code      
import syscadif
import numpy as np
import matplotlib.pyplot as plt
 
sc = syscadif.SysCADCom()
ScdDir = r'C:\SysCAD%d' % sc.build
ScdPrj = r'\Examples\65 Smelting\Copper Flash Furnace Example.spf\Project.spj'
sc.OpenProject(ScdDir+ScdPrj)  ## Open project 

VarTags = {"Cu_Concentrate.QmReqd (t/h)": (62, 100, 2), 
           "Controllers.Cfg.[2].Spt":     (1240, 1320, 4),
           ## add further variables as needed Tagstr: (vmin, vmax, npts)
           }
SensTags = ["P_005.Qm (t/h)", "$Solver.Convergence.Iterations"]

V  = tuple(np.linspace(*v) for v in VarTags.values())
mgrid  = np.meshgrid(*V)

@np.vectorize
def f(*U):
    '''Funtion to set primary variables, run SysCAD, and return sensitivity variable'''
    for u, t in zip(U, VarTags.keys()):
        sc[t] = u
    sc.run()
    return tuple(sc[t] for t in SensTags)
 
## Now evaluate the SysCAD model over a n-dimensional grid (n is the number of VarTags)
## This returns a tuple of n-dimensional arrays, each being the values of the sensitivity variables at the gridpoints defined by meshgrid.
Z = f(*mgrid)  

X, Y = mgrid
fig = plt.figure()
ax = fig.add_subplot(1,2,1, projection='3d')
ax.plot_wireframe(X, Y, Z[0])
ax.set_xlabel("Primary Var")
ax.set_ylabel("Secondary Var")
ax.set_zlabel("Sensitivity Var")
ax = fig.add_subplot(1,2,2, projection='3d') 
ax.plot_wireframe(X, Y, Z[1])
ax.set_xlabel("Concentrate Flow (tph)")
ax.set_ylabel("Furnace Temperature (C)")
ax.set_zlabel("Iterations")
 
plt.show()
sc.Close()
CuSensitivity.png CuSensitivity1.png

If we just want to create a 2D plot like the original example (see right image above)

for Y, v in zip(Z[0], V[1]):
    plt.plot(V[0], Y, label = ("%6.1f" % v))
plt.legend()
plt.title("Sensitivity Analysis")
plt.ylabel("Sensitive Variable")
plt.xlabel("Primary Variable")

This example can easily be extended to more primary tags and sensitivity variables. Just add more items to the dictionary of primary tags (VarTags), and list of sensitivity tags. The code creates a meshgrid of the ranges, and call the vectorized function f on the meshgrid. This returns a tuple of the sensitivity variables as numpy arrays, which can be used for plotting or further analysis. This may involve a lot of evaluations of SysCAD!

If we had five primary variables, each to be evaluated at four points and returning 10 sensitivity variables, we would require [math]\displaystyle{ 4^5 = 1024 }[/math] evaluations of the flowsheet, and the vectorized function call would return 10240 data points.

Plotting PSD Data

Plotting Data in Log Scale

In this next example, we’ve copied seven columns of data. Column 1 represents the size fractions. Columns 2 to 4 contain the particle size distributions (PSD) for the Feed, Overflow (OF), and Underflow (UF) streams. Columns 5 to 7 present the Fe₂O₃ PSD for the Feed, OF, and UF respectively.

import numpy as np
from io import StringIO
import win32clipboard as wcb
import matplotlib.pyplot as plt

# Read data from clipboard
wcb.OpenClipboard()
data = wcb.GetClipboardData()
wcb.CloseClipboard()

# Convert clipboard string data to a NumPy array using StringIO
data_io = StringIO(data)
lines = data_io.readlines()

# Extract header and data
header = lines[0].strip().split('\t')  # Adjust delimiter if needed
data_array = np.genfromtxt(StringIO(''.join(lines[1:])), delimiter='\t')

# Extract x values (first column) and y values (remaining columns)
x = data_array[:, 0]
y_data = data_array[:, 1:]

# Plot each data series with legend from header
plt.figure(figsize=(10, 6))
for i in range(y_data.shape[1]):
    if i >= 3: # Series 4, 5, 6 (0-based index)
        plt.plot(x, y_data[:, i], label=header[i+1], marker='o') # round bullet marker
    else:
        plt.plot(x, y_data[:, i], label=header[i+1]) # regular line

# Set log scale for x-axis
plt.xscale('log')

# Label axes and add title
plt.xlabel('Size (mm)')
plt.ylabel('Cumulative Fraction Passing (%)')
plt.title('Log Scale Plot of Cumulative Fraction Passing')

# Add legend and grid
plt.legend()
plt.grid(True, which="both", linestyle="-", color='lightgrey')

# Display the plot
plt.tight_layout()
plt.show()
Pypsd2.png

Sub Plots

Size distribution data can vary widely, and when dealing with multiple ore types and process streams, there are many ways to visualize it. In this example, we have three columns of data—each representing size measurements for an ore species across three streams. The dataset contains 60 rows in total, with 20 data points per stream. We can visualize this data either by stream or by ore species.

import numpy as np
from io import StringIO
import win32clipboard as wcb
import matplotlib.pyplot as plt
wcb.OpenClipboard()  ## data in clipboard as columns for each stream
data = wcb.GetClipboardData()
wcb.CloseClipboard()

X, Y, Z = np.loadtxt(StringIO(data), skiprows = 1, unpack=True) 
X = X.reshape((3, -1))
Y = Y.reshape((3, -1))
Z = Z.reshape((3, -1))

plt.ion()
fig, ax = plt.subplots(3, 2, sharex = True, sharey = True)
titles = ["Fe2O3(s)","Fe3O4(s)","SiO2(s)",]
cols = ["r-", 'g-','b-'] 
for i, t in zip(range(3), titles):#
  for a, c in zip([X, Y, Z], ["r-", 'g-','b-']):
    ax[i, 0].plot(a[i], c)#
  ax[i, 0].set_title(t)

titles = ["Feed","Underflow","Overflow",] 
for i, a, t in zip(range(3), [X, Y, Z], titles):
  for j in range(3):
    ax[i, 1].plot(a[j])	
  ax[i, 1].set_title(t)

plt.tight_layout()
plt.show()
Pypsd.png

Plotting in 3D

This example pulls particle size distributions data for a series of pipes in the Distributed Precipitation and Classification Example without ever running SysCAD. It then makes a 3d plot showing the variation in PSD curve along the tank row.

import syscadif  ## SysCAD COM interface class
import numpy as np
import matplotlib.pyplot as plt

sc = syscadif.SysCADCom()
ScdDir = r'C:\SysCAD%d' % sc.build
ScdPrj = r'\Examples\10 Alumina\PSD Precipitation & Classification Example.spf\Project.spj'
 
sc.OpenProject(ScdDir+ScdPrj)  ## Open project 
pipelist = ["P_115_", "P_116_"]+["P_0%02d" % i for i in range(4, 24, 2)]+["P_121"]
sizetags = [".Qo.Sz.I%d.FP.Al[OH]3(s) (%%)" % i for i in range(28)]

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

X, Y = np.meshgrid(np.arange(len(pipelist)),
                  np.arange(len(sizetags)))
Z = np.array([ [sc[pipe+s] for pipe in pipelist] for s in sizetags])

ax.plot_wireframe(X, Y, Z)
plt.show()
Comexample.png

Plotting Thermophysical Data

We can generate thermophysical data using Species Property reports when viewing species data: here we select all the different phases of a particular species and plot the specific heat, entropy, enthalpy and Gibbs free energy. The data is put in the clipboard, and it is easy to access the data directly from the clipboard and plot it:

This script demonstrates accessing the clipboard, using numpy to read blocks of text data into arrays, and then plotting this, along with labels and legends.

import win32clipboard as wcb
import numpy as np
import matplotlib.pyplot as plt
from io import StringIO
import os

blanklinesep=os.linesep+os.linesep    ## Empty line for split
  
def readCB():
    #'''Read data from the clipboard'''
    wcb.OpenClipboard()
    data = wcb.GetClipboardData()
    wcb.CloseClipboard()
    return data.strip().split(blanklinesep)  ## split on empty lines  
 
def main():
    #'''Each data block has three header lines indicating property, pressure and species'''
    pd = readCB()
    for i, z  in enumerate(pd):
        s = StringIO(z)     ## StringIO makes a 'file like' object from a string
        yvar = s.readline().strip()  ## The independent variable
        pvar = s.readline()  ## the pressure
        vvars = s.readline().strip().split("\t")   ## Temperature units and species
        ddata = np.loadtxt(s, unpack = True)
        xdata, ydata = ddata[0], ddata[1:]
        plt.figure(i+1)
        for j, yy in enumerate(ydata):
            plt.plot(xdata, yy, label = vvars[j+1])
            plt.xlabel(vvars[0])
            plt.ylabel(yvar)
            plt.legend()
    plt.show()
 
if __name__ == "__main__":
    main()
Fe4.png

In particular the Gibbs function curves for the phases intersect at the melting and boiling points. At low temperatures the solid (cr) phase has the lowest free energy, and as the temperature increases, the liquid and then vapour phases have lower free energy.

Python10.png Pbfe.png

Plotting Species Flows

As well as the usual plots we can create bar or pie charts. The Qo tab on any stream has the option to copy the species to the clipboard, so we can create a chart showing the species flows.

from io import StringIO
import matplotlib.pyplot as plt
import numpy as np
import win32clipboard as wcb
 
wcb.OpenClipboard()
cbtxt = wcb.GetClipboardData()
wcb.CloseClipboard()

res = np.loadtxt(StringIO(cbtxt), dtype = {'names': ('species', 'unit', 'value'),
                                           'formats' : ('U30', 'U10', 'f8')}, skiprows = 2)
v = res['value']
sp = res['species'][v>0]   ## species with non-zero flows
fig, ax = plt.subplots()
y_pos = np.arange(len(sp))
ax.set_xscale('log')
ax.barh(y_pos, v[v>0], align='center')
ax.set_yticks(y_pos)
ax.set_yticklabels(sp)
ax.invert_yaxis()  # labels read top-to-bottom
ax.set_xlabel('Mass Flow')
ax.set_title('Species plot')
 
plt.show()
Pyspecies.png
  • This shows using structured arrays in numpy, which is useful when we have mixed data types, in this case string data for the species, and numerical data for the flows.
  • The notation v[v>0] is showing use of an index array.
  • v>0 is an array of flags indicating which elements of the array v are non-zero.
>>> v>0
array([ True,  True,  True,  True,  True,  True,  True,  True, False,
       True, False,  True, False, False, False, False, False, False,
      False, False, False, False, False, False, False, False, False,
      False, False, False, False, False, False,  True, False, False,
      False, False, False, False, False, False, False])
  • When we use an array as an index, it returns just the elements for which the index array are True, i.e. the non-zero elements of the original array.
 >>> v[v>0]
 array([3.75694246e+02, 1.76991503e+02, 2.72699554e+01, 1.28469206e+01,
      6.32469251e-01, 2.60232738e+00, 5.87385513e-01, 1.45618792e+01,
      5.32841733e-02, 2.42947148e+01, 4.20754018e-02])
  • Similarly res['species'] is an array of strings having the species names, while res['species'][v>0] is an array of names of species with non-zero flow.
  • We plot this on a logarithmic scale to show relative magnitudes of the flows of the species.