FMU - Approximation Influence



In [1]:
# Import the needed Packages
# FMU Simulation
import MoBASimulator as mb
# Numpy
import numpy as np
# Bokeh for Plotting
#import bokeh.plotting as bk
#import bokeh.io as bi
#from bokeh.io import export_svgs
#bi.output_notebook()
# select a palette
#from bokeh.palettes import Spectral10  as palette
#from bokeh.models import Legend, LabelSet, Label, BoxAnnotation
# Matplotlib for Plotting
import matplotlib.pyplot as plt

# Algorithms
import Algorithms as alg
from Algorithms import TUBScolorscale
from Algorithms import cm2in
# Make the relevant Inputs
export_folder = "../Latex/Graphics/"
model_name = "Approximation_Limits"

# Plot width and heigth in cm
plot_width = 19.
plot_height = 7.5

# Model Parameter
k = [[1,.3],[.3,1]]
t = [[99.,101.],[102.,98.]]
l = [[.1,.1],[50,.1]]
# Rosenbrock
#k = [[1.,2./3.],[1.,1.]]
#t = [[1.,1./3.],[1.,1.]]
#l = [[1e-10,1e-10],[1e-10,1e-10]]

# Woodberry
#k = [[12.8,-18.9],[6.6,-19.4]]
#t = [[16.7,21.0],[10.9,14.4]]
#l = [[1,3],[7,3]]



# Make variable storage
yl = []
ul = []
timel = []
TApprox = np.zeros((2,2,10))
LApprox = np.zeros((2,2,10))
TReal = np.zeros((2,2,10))
LReal = np.zeros((2,2,10))
# Show log window
sim = mb.Simulator()
sim.showLogWindow()

for delay in range(10,105,5):
    # Make variable storage
    yex = []
    uex = []
    timex = []
    
    # The needed Parameter
    K = np.zeros((2,2))
    T = np.zeros((2,2))
    L = np.zeros((2,2))

    # Load a Model
    sim.clear()
    sim.loadModel("C:/Users/juliu/Documents/Thesis/Modelica/FMU/2_2/Masterthesis_Models_mimo_0processmodel.fmu")

    # Parameter Values
    params = {}
    # Loop over system size
    for outputs in range(1,3):
        for inputs in range(1,3):
            # Process Parameter 
            # System Gain
            params.update({"fmu.num["+str(outputs)+","+str(inputs)+",1]": k[outputs-1][inputs-1]})
            # System Lag
            params.update({"fmu.den["+str(outputs)+","+str(inputs)+",2]": 1})
            # System Lag
            params.update({"fmu.den["+str(outputs)+","+str(inputs)+",2]": 1})
            # System Delay
            params.update({"fmu.delay["+str(outputs)+","+str(inputs)+"]": (delay/100.0)*l[outputs-1][inputs-1]})
            params.update({"fmu.den["+str(outputs)+","+str(inputs)+",1]":t[outputs-1][inputs-1] - (delay/100.0)*l[outputs-1][inputs-1]})

    # Set Parameter and show for checking
    sim.set(params)
    #sim.showParameterDialog()
    #print(params)



    # First run, Input 1 -> Output 1 & 2
    sim.set({"fmu.u[1]": 1,"fmu.u[2]": 0})
    # Set timestep = 1e-2, endtime = 100
    sim.resetModelState()
    res = sim.simulate(0.01, 500)

    # Get the signals
    y = res["fmu.y[1]"]
    y2 = res["fmu.y[2]"]
    u = res["fmu.u[1]"]
    time = res["time"]

    # Get TF from Input 1 to Output 1
    K[0][0],T[0][0],L[0][0]=alg.Integral_Identification(y,u,time)
    # Get TF from Input 1 to Output 2
    K[1][0],T[1][0],L[1][0]=alg.Integral_Identification(y2,u,time)


    # Second run, Input 2 -> Output 1 & 2
    # Reload the Model to set everything to zero
    sim.resetModelState()
    sim.set({"fmu.u[1]":0, "fmu.u[2]":1})
    # Set timestep = 1e-2, endtime = 100
    res=sim.simulate(0.01, 500)

    # Get the signals
    y = res["fmu.y[1]"]
    y2 = res["fmu.y[2]"]
    u = res["fmu.u[2]"]
    time = res["time"]

    # Get TF from Input 2 to Output 1
    K[0][1],T[0][1],L[0][1] = alg.Integral_Identification(y,u,time)
    # Get TF from Input 2 to Output 2
    K[1][1],T[1][1],L[1][1] = alg.Integral_Identification(y2,u,time)


    # Loop over the methods for controller design
    for methods in range(0,3):
        if methods == 0:
            # Make a decentralized controller
            KY,B,D = alg.Control_Decentral(K,T,L)
        elif methods == 1:
            # Make a decoupled controller with Astrom
            KY,B,D = alg.Control_Astrom(K,T,L,np.eye(2,2))
        else:
            # Make a decoupled controller with modified Astrom
            KY,B,D = alg.Control_Decoupled(K,T,L,np.eye(2,2))

        # Make zeros to 1e-10 for numerical stable process
        KY[KY==0] = 1e-10
        B[B==0] = 1e-10
        D[D==0] = 1e-10


        # Load the closed loop model
        sim.clear()
        sim.loadModel("C:/Users/juliu/Documents/Thesis/Modelica/FMU/2_2/Masterthesis_Models_mimo_0closedloop.fmu")
        # Parameter Valuess
        params = {}

        # Loop over system size
        for outputs in range(1,3):
            for inputs in range(1,3):
                # Controller Parameter 
                # Proportional Gain
                params.update({"fmu.kp["+str(outputs)+","+str(inputs)+"]": KY[outputs-1][inputs-1][0]})
                 # Integral Gain
                params.update({"fmu.ki["+str(outputs)+","+str(inputs)+"]": KY[outputs-1][inputs-1][1]})
                # Decoupler
                params.update({"fmu.d["+str(outputs)+","+str(inputs)+"]": D[outputs-1][inputs-1]})
                # Set Point Weight
                params.update({"fmu.b["+str(outputs)+","+str(inputs)+"]": B[outputs-1][inputs-1]})
                # Process Parameter 
                # Process Parameter 
                # System Gain
                params.update({"fmu.num["+str(outputs)+","+str(inputs)+",1]": k[outputs-1][inputs-1]})
                # System Lag
                params.update({"fmu.den["+str(outputs)+","+str(inputs)+",2]": 1})
                # System Delay
                params.update({"fmu.delay["+str(outputs)+","+str(inputs)+"]": (delay/100.0)*l[outputs-1][inputs-1]})
                params.update({"fmu.den["+str(outputs)+","+str(inputs)+",1]":t[outputs-1][inputs-1] - (delay/100.0)*l[outputs-1][inputs-1]})
            
        # Set Parameter Values
        sim.set(params)

        # First run, Input 1 -> Output 1 & 2
        sim.resetModelState()
        sim.set({"fmu.u[1]": 1.,"fmu.u[2]": 0.})
        res=sim.simulate(1, 1500)
        # Second run, Input 2-> Output 1 & 2
        sim.set({"fmu.u[1]": 1.,"fmu.u[2]": 1})
        res=sim.simulate(1, 3000)

        # Store the output
        yex.append([res["fmu.y[1]"],res["fmu.y[2]"]])
        uex.append([res["fmu.u[1]"],res["fmu.u[2]"]])
        timex.append([res["time"],res["time"]])
        
    yl.append(yex)
    ul.append(uex)
    timel.append(timex)


Aström Detuning Iterations:0
Aström Detuning Iterations:0
Modified Detuning Iterationts 0
Modified Detuning Iterationts 0
Aström Detuning Iterations:0
Aström Detuning Iterations:0
Modified Detuning Iterationts 0
Modified Detuning Iterationts 0
Aström Detuning Iterations:0
Aström Detuning Iterations:0
Modified Detuning Iterationts 0
Modified Detuning Iterationts 0
Aström Detuning Iterations:0
Aström Detuning Iterations:0
Modified Detuning Iterationts 0
Modified Detuning Iterationts 0
Aström Detuning Iterations:0
Aström Detuning Iterations:0
Modified Detuning Iterationts 0
Modified Detuning Iterationts 0
Aström Detuning Iterations:0
Aström Detuning Iterations:0
Modified Detuning Iterationts 0
Modified Detuning Iterationts 0
Aström Detuning Iterations:0
Aström Detuning Iterations:0
Modified Detuning Iterationts 0
Modified Detuning Iterationts 0
Aström Detuning Iterations:0
Aström Detuning Iterations:0
Modified Detuning Iterationts 0
Modified Detuning Iterationts 0
Aström Detuning Iterations:0
Aström Detuning Iterations:0
Modified Detuning Iterationts 0
Modified Detuning Iterationts 0
Aström Detuning Iterations:0
Aström Detuning Iterations:0
Modified Detuning Iterationts 0
Modified Detuning Iterationts 0
Aström Detuning Iterations:0
Aström Detuning Iterations:0
Modified Detuning Iterationts 0
Modified Detuning Iterationts 0
Aström Detuning Iterations:0
Aström Detuning Iterations:0
Modified Detuning Iterationts 0
Modified Detuning Iterationts 0
Aström Detuning Iterations:0
Aström Detuning Iterations:0
Modified Detuning Iterationts 0
Modified Detuning Iterationts 0
Aström Detuning Iterations:0
Aström Detuning Iterations:0
Modified Detuning Iterationts 0
Modified Detuning Iterationts 0
Aström Detuning Iterations:0
Aström Detuning Iterations:0
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-1-3fe4644df1a4> in <module>()
    188 
    189         # Store the output
--> 190         yex.append([res["fmu.y[1]"],res["fmu.y[2]"]])
    191         uex.append([res["fmu.u[1]"],res["fmu.u[2]"]])
    192         timex.append([res["time"],res["time"]])

KeyError: 'fmu.y[1]'

In [ ]:
plt.rcParams['svg.fonttype'] = 'none'
plt.style.use('seaborn-whitegrid')

fig, ax = plt.subplots(2, sharex=True,figsize = cm2in(plot_width,2.*plot_height))
for i in range(0,10):
#ax[0].plot(res["time"],res["fmu.y[1]"], color = TUBScolorscale[0], linewidth = 1., linestyle = "dashed")
    time = timel[i]
    y = yl[i]
    ax[0].plot(time[0][0],y[0][0], color = TUBScolorscale[9], linewidth = 1.)
    ax[0].plot(time[1][0],y[1][0], color = TUBScolorscale[1], linewidth = 1.)
    ax[0].plot(time[2][0],y[2][0], color = TUBScolorscale[3], linewidth = 1.)
    ax[0].grid(True)
    ax[0].set_ylabel('$y_1$')

    #ax[1].plot(res["time"],res["fmu.y[2]"], color = 'k', linewidth = 1., linestyle = "dashed", label="Reference")
    ax[1].plot(time[0][1],y[0][1], color = TUBScolorscale[9], linewidth = 1.)
    ax[1].plot(time[1][1],y[1][1], color = 'grey', linewidth = 1., alpha = 0.6)
    ax[1].plot(time[2][1],y[2][1], color = TUBScolorscale[3], linewidth = 1.)
    ax[1].grid(True)
    ax[1].set_ylabel('$y_2$')
    ax[1].set_xlabel('Time [s]')

time = timel[0]
y = yl[0]
ax[1].plot(time[1][1],y[1][1], color = TUBScolorscale[5], linewidth = 1., linestyle="dashed", label = 'L_{22} = 5')

time = timel[9]
y = yl[9]
ax[1].plot(time[1][1],y[1][1], color = TUBScolorscale[1], linewidth = 1., linestyle="dashed", label = 'L_{22} = 50')


plt.legend(loc="lower right")
plt.savefig(export_folder+model_name+".svg")
plt.show()

In [ ]:
# Bokeh for Plotting
import bokeh.plotting as bk
import bokeh.io as bi
from bokeh.layouts import gridplot

p1 = bk.figure(title = "Influence of the Approximation", height= 300, width = 800)
p2 = bk.figure(height= 300, width = 800)
for i in range(0,10):
    y = yl[i]
    time = timel[i]
    p1.line(time[0][0],y[0][0], color = "red")
    p1.line(time[1][0],y[1][0], color = "blue")
    p1.line(time[2][0],y[2][0], color = "green")

    p2.line(time[0][1],y[0][1], color = "red")
    p2.line(time[1][1],y[1][1], color = "blue")
    p2.line(time[2][1],y[2][1], color = "green")
p = gridplot([[p1],[p2]])
bk.show(p)

In [ ]: