Commit 8a909991 authored by Blais, Chris's avatar Blais, Chris
Browse files

adding UCI operating point notebook and other materials

parent e35f0ea5

documenting_failures/get_nonlinear_integer_x0_membrane.ipynb

0 → 100644
+323 −0
Original line number Diff line number Diff line 
%% Cell type:code id:1c623c5c tags:



``` python

import sympy as sy

import numpy as np

import os

import sys

import scipy as sp



from itertools import product



from sympy.core.numbers import int_valued



module_path = os.path.join(os.path.dirname(os.path.abspath('')))

if module_path not in sys.path:

    sys.path.insert(0, module_path)

else:

    print("path already in sys.path")



import matplotlib.pyplot as plt

from numerical_labelling.labelling import *

```



%% Cell type:markdown id:20ad12df tags:



### simplified membrane example:

a simple layout (1 membrane), 3*3=9 measurements



- 1 hydraulic balance

- 1 component balance

- 1 pressure constraint

- 1 concentration constraint



%% Cell type:code id:5c7d0989 tags:



``` python

def add_sltn(sltn_dict, varb_name, constr):

    """add a solved equation to sltn dict

    example:

    varb_name = x3

    constr = m1x1 + m2x2 - m3x3

    sltn_dict.update({'x3': (m1x1 + m2x2)/m3})

    """

    sltn = sy.solve(constr, varb_name)

    if len(sltn) !=1:

        raise ValueError(f"multiple solutions when trying to solve for {varb_name}")



    if varb_name in sltn_dict.keys():

        print(f"variable {varb_name} in eq {constr} already solved for in eq {sltn_dict[varb_name]}")

        # raise KeyError(f"{varb_name} already solved for")

    sltn_dict.update({varb_name:sltn[0]})



    return sltn_dict

```



%% Cell type:code id:83011f78 tags:



``` python

def add_constraint(constr_dict, eqn):

    """

    add constraint to constraint dict

    """

    indx = len(constr_dict)

    constr_dict[indx] = eqn

    return constr_dict



def constr_split(

    min, mout1, mout2,

    xin, xout1, xout2,

    pin, pout1, pout2,

    constr, sltns=None,

    ):

    """

    constraint for splitter unit

    """

    constr_hydraulic = min - mout1 - mout2

    constr_eq_conc1 = xin - xout1

    constr_eq_conc2 = xin - xout2

    constr_press1 = pin - pout1

    constr_press2 = pin - pout2

    constr = add_constraint(constr, constr_hydraulic)

    constr = add_constraint(constr, constr_eq_conc1)

    constr = add_constraint(constr, constr_eq_conc2)

    constr = add_constraint(constr, constr_press1)

    constr = add_constraint(constr, constr_press2)

    if sltns is not None:

        sltns = add_sltn(sltns, xout1.name, constr_eq_conc1)

        sltns = add_sltn(sltns, xout2.name, constr_eq_conc2)

        sltns = add_sltn(sltns, pout1.name, constr_press1)

        sltns = add_sltn(sltns, pout2.name, constr_press2)

        return constr, sltns

    return constr



def constr_mix(

        min1, min2, mout,

        xin1, xin2, xout,

        pin1, pin2, pout,

        constr, sltns=None

        ):

    """

    constraint for 2 pipes coming together

    """

    constr_hydraulic = min1 + min2 - mout

    constr_component = min1*xin1 + min2*xin2 - mout*xout

    constr_press1 = pin1 - pout

    constr_press2 = pin2 - pout

    constr = add_constraint(constr, constr_hydraulic)

    constr = add_constraint(constr, constr_component)

    constr = add_constraint(constr, constr_press1)

    constr = add_constraint(constr, constr_press2)



    if sltns is not None:

        sltns = add_sltn(sltns, xout.name, constr_component)

        sltns = add_sltn(sltns, pout.name, constr_press1)

        sltns = add_sltn(sltns, pin2.name, constr_press2)

        return constr, sltns

    return constr

```



%% Cell type:code id:bbd31430 tags:



``` python

def constr_mex(

        mfeed, mperm, mret,

        xfeed,xperm, xret,

        pfeed, pperm, pret,

        A, B, C, D,

        constr, sltns = None

    ):

    # # see if multiple a, b, c, and ds defined

    # A, B, C, D = sy.symbols("A, B, C, D")

    # varb_symbs = [A, B, C, D]

    # labels = [ "A", "B", "C", "D"]

    # varbs.update(dict(zip(labels, varb_symbs)))



    constr_hyd = mfeed - mperm  - mret

    constr_component = xfeed*mfeed -xperm*mperm -xret*mret

    constr_osm_press = mperm - A*((pfeed - pperm) - C*(xfeed-xperm))

    constr_diff_conc = mperm*xperm - B*C*(xfeed-xperm) # xperm may be very small, but for integer example maybe it's better to keep

    constr_major_loss = pfeed - pret - D*(mfeed+mret)**2 # major loss from feed to concentrate side.

    constr = add_constraint(constr, constr_hyd)

    constr = add_constraint(constr, constr_component)

    constr = add_constraint(constr, constr_osm_press)

    constr = add_constraint(constr, constr_diff_conc)

    constr = add_constraint(constr, constr_major_loss)



    # return solutions for the membrane constraint equations

    if sltns is not None:

        sltns = add_sltn(sltns, xperm.name, constr_diff_conc)

        sltns = add_sltn(sltns, pperm.name, constr_osm_press)

        sltns = add_sltn(sltns, pret.name, constr_major_loss)

        return constr, sltns



    return constr

```



%% Cell type:code id:fe297a26 tags:



``` python

def constr_pump(

        min, mout,

        xin, xout,

        pin, pout,

        E, F, G,

        constr, sltns=None,

    ):

    """

    E: second order coefficient (Ex^2)

    F: first order coefficient (Fx)

    G: intercept for pump curve

    """

    constr_mass = mout - min

    constr_conc = xout - xin

    constr_pump = (pout - pin) - (E*min**2 + F*min + G)

    constr = add_constraint(constr, constr_mass)

    constr = add_constraint(constr, constr_conc)

    constr = add_constraint(constr, constr_pump)



    if sltns is not None:

        sltns = add_sltn(sltns, xout.name, constr_conc)

        sltns = add_sltn(sltns, pout.name, constr_pump)

        return constr, sltns

    return constr



```



%% Cell type:code id:1c04882e tags:



``` python

def constr_valve(

        min, mout,

        xin, xout,

        pin, pout,

        K,

        constr, sltns= None,

    ):

    # # see if multiple a, b, c, and ds defined

    # A, B, C, D = sy.symbols("A, B, C, D")

    # varb_symbs = [A, B, C, D]

    # labels = [ "A", "B", "C", "D"]

    # varbs.update(dict(zip(labels, varb_symbs)))



    constr_mass = mout - min

    constr_conc = xout - xin

    constr_major_loss = pin - pout - K*(min+mout)**2 # major loss from feed to concentrate side.



    constr = add_constraint(constr, constr_mass)

    constr = add_constraint(constr, constr_conc)

    constr = add_constraint(constr, constr_major_loss)



    if sltns is not None:

        sltns = add_sltn(sltns, xout.name, constr_conc)

        sltns = add_sltn(sltns, pout.name, constr_major_loss)

        return constr, sltns

    return constr

```



%% Cell type:code id:1c261e9a tags:



``` python

constr_ex = {}

m1, m2, m3 = sy.symbols("m1 m2 m3")

x1, x2, x3 = sy.symbols("x1 x2 x3")

p1, p2, p3 = sy.symbols("p1 p2 p3")

A1, B1, C1, D1 = sy.symbols("A1 B1 C1 D1")

sy_varbs = [

    m1, m2, m3,

    x1, x2, x3,

    p1, p2, p3,

    A1, B1, C1, D1,

]  # List of symbolic variables

sltns_ex = {}

varbs_ex = {varb.name:varb for varb in sy_varbs}

constr_ex, sltns_ex = constr_mex(

    mfeed=m1, mperm=m2, mret=m3,

    xfeed=x1, xperm=x2, xret=x3,

    pfeed = p1, pperm=p2, pret=p3,

    A=A1, B=B1, C=C1, D=D1,

    constr=constr_ex, sltns=sltns_ex,

)

```



%% Cell type:code id:4a211fcc tags:



``` python

sltns_ex

```



%% Output



    {'x2': B1*C1*x1/(B1*C1 + m2),

     'p2': -C1*x1 + C1*x2 + p1 - m2/A1,

     'p3': -D1*(m1 + m3)**2 + p1}



%% Cell type:code id:c977f98e tags:



``` python

# solve membrane concentrate passage for x2

x2_constr = sy.solve(constr_ex[3], varbs_ex["x2"])[0]

print("membrane dp constr")

display(constr_ex[2])

print("membrane pipe fitting loss term")

display(constr_ex[4])

p1_sub1 = sy.solve(constr_ex[2], varbs_ex['p2'])[0]

p1_sub2 = sy.solve(constr_ex[4], varbs_ex['p3'])[0]

display(p1_sub1)

display(p1_sub2)

# p1_constr = p1_sub1 - p1_sub2

# display(p1_constr)

# p3_constr = sy.solve(p1_constr, varbs_ex['p3'])[0]

# display(p3_constr)

```



%% Output



    membrane dp constr



    $\displaystyle - A_{1} \left(- C_{1} \left(x_{1} - x_{2}\right) + p_{1} - p_{2}\right) + m_{2}$



    membrane pipe fitting loss term



    $\displaystyle - D_{1} \left(m_{1} + m_{3}\right)^{2} + p_{1} - p_{3}$



    $\displaystyle - C_{1} x_{1} + C_{1} x_{2} + p_{1} - \frac{m_{2}}{A_{1}}$



    $\displaystyle - D_{1} \left(m_{1} + m_{3}\right)^{2} + p_{1}$



%% Cell type:code id:c2d652d3 tags:



``` python

x0 = dict(zip(list(varbs_ex.keys()), [1] * len(varbs_ex)))

x0['m1'] = 100

x0['m2'] = 20

x0['m3'] = (x0['m1'] - x0['m2'])

x0['x1'] = 30*7

x0['x2'] = x2_constr.subs(x0)

x0['x3'] = ((x0['m1']*x0['x1'] - x0['m2']*x0['x2'])/x0['m3'])

x0['p1'] = 80 + 32340

x0['p2'] = p1_sub1.subs(x0)

x0['p3'] = p1_sub2.subs(x0)

x0['A1'] = -1

x0['B1'] = 25

x0['C1'] = 1

x0['D1'] = 1

x0



```



%% Output



    {'m1': 100,

     'm2': 20,

     'm3': 80,

     'x1': 210,

     'x2': 10,

     'x3': 260,

     'p1': 32420,

     'p2': 32200,

     'p3': 20,

     'A1': -1,

     'B1': 25,

     'C1': 1,

     'D1': 1}

documenting_failures/get_nonlinear_integer_x0_pilot.ipynb

0 → 100644
+1289 −0

File added.

Preview size limit exceeded, changes collapsed.

fig/pilot_layout.png

0 → 100644
+129 KiB
Loading image diff...

numerical_labelling/plot_layout.py

+77 −0
Original line number Diff line number Diff line
@@ -42,6 +42,83 @@ def plot_bl_layout(layout_plt, varbs_bl): 
                ) # m1

            



def plot_layout_any(layout_plt, layout_name, varbs, x_coords = [], y_coords = [], ):



    """

    layout_name: the name of the layout, 'pilot', 'minimal'

    """

    img_folder = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), "fig")

    # print(img_folder)

    mem_img_path = os.path.join(img_folder, f"{layout_name}_layout.png")

    mem_img = np.asarray(Image.open(mem_img_path))

    fig, ax = plt.subplots(figsize=(45,20))

    ax.imshow(mem_img, aspect="auto")



    conc_xoffset = 120 # x offset for concentration

    conc_yoffset = 0

    press_xoffset = 220 # offset for pressure

    press_yoffset = 0





    # xy positions for text

    match layout_name:

        case "example":

            xs = [200, 700, 1050, 1800, 1575, 1850, 1300, 950]

            ys = [200, 200, 200,   200,  850,  800,  800, 850]



        case "pilot":

            if len(x_coords) > 0:

                xs = x_coords

            else:

                xs = [

                    250, 750, 1200, 1425, 1700, 2250, 2350, 2250, 2800, 2900, 2800, 3350, 3550, 3350, 3550, 1800, 650

                ]

            if len(y_coords):

                ys = y_coords

            else:

                ys =  [

                    450, 450, 450, 200, 1075, 500, 175, 1075, 500, 175, 1075, 500, 175, 1250, 1475, 1475, 1475

                ] 

    

    if len(varbs) == len(xs):

        mcp = True, False, False

    elif len(varbs) == 2*len(xs):

        mcp = True, True, False

    elif len(varbs) == 3*len(xs):

        mcp = True, True, True

    

    # append coordinates for x and y positions for conc, press

    xs_all = copy.deepcopy(xs)

    ys_all = copy.deepcopy(ys)

    if mcp[1]:

        for x in xs: 

            xs_all.append(x+conc_xoffset)

        for y in ys:

            ys_all.append(y+conc_yoffset)

    if mcp[2]:

        for x in xs: 

            xs_all.append(x+press_xoffset)

        for y in ys:

            ys_all.append(y+press_yoffset)





    names = list(varbs.keys())



    for i, coords in enumerate(zip(xs_all, ys_all)):

        if layout_plt[i]:

            if names[i].lower().startswith("m"):

                color = "blue"

            elif names[i].lower().startswith("x"):

                color = "green"

            elif names[i].lower().startswith("p"):

                color = "red"



            ax.text(coords[0], coords[1], names[i], fontsize=30, weight="bold", color = color) # m1

        else: 

            ax.text(

                coords[0], coords[1], names[i], fontsize=30, color = "grey", alpha = 0.5,

                ) # m1



def make_pretty_labels(varbs, layout):



    """