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()
# Algorithms
import Algorithms_Graphic as alg

In [2]:
# Model Parameter
k = [[10.,.5],[5.,8.]]
l = [[1e-10,1e-10],[1e-10,1e-10]]
# Generate stable polynomial with roots between -1/1....-1/20 rad/s 
r = np.random.uniform(-1.0,-8.,(2,2,8))
t = np.zeros((2,2,9))
for outputs in range(0,2):
    for inputs in range(0,2):
        t[outputs,inputs,:] = np.polynomial.polynomial.polyfromroots(r[outputs,inputs,:])


# The needed Parameter
K = np.zeros((2,2))
T = np.zeros((2,2))
L = np.zeros((2,2))

In [3]:
# Load a Model
sim = mb.Simulator()
sim.clear()
sim.loadModel("C:/Users/juliu/Documents/Thesis/Modelica/FMU/2_2_n9/Masterthesis_Models_mimo_0processmodel.fmu")
sim.setOperationMode('FMU for ModelExchange')

# Show log window
sim.showLogWindow()

# 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
        for order in range(0,9):
            params.update({"fmu.den["+str(outputs)+","+str(inputs)+","+str(order+1)+"]": t[outputs-1][inputs-1][order]})
        #params.update({"fmu.den["+str(outputs)+","+str(inputs)+",9]": 1})
        # System Delay
        params.update({"fmu.delay["+str(outputs)+","+str(inputs)+"]": l[outputs-1][inputs-1]})

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

In [4]:
# First run, Input 1 -> Output 1 & 2
sim.set({"fmu.u[1]": 1,"fmu.u[2]": 0})
# 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[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, graphics='on')
# Get TF from Input 1 to Output 2
K[1][0],T[1][0],L[1][0]=alg.Integral_Identification(y2,u,time, graphics='on')


# 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, graphics='on')
# Get TF from Input 2 to Output 2
K[1][1],T[1][1],L[1][1] = alg.Integral_Identification(y2,u,time, graphics='on')
K,T,L


Out[4]:
(array([[ 10.00000154,   0.5000084 ],
        [  5.00002465,   8.00006863]]), array([[ 17.42971964,  17.36779594],
        [ 53.70478135,  13.94935342]]), array([[ 27.66291934,  26.86354554],
        [  1.19638372,  21.35359732]]))

In [5]:
# Make a controller
KY,B,D = alg.Control_Decentral(K,T,L)
#KY,B,D = alg.Control_Astrom(K,T,L,0.01*np.eye(2,2))
#KY,B,D = alg.Control_Decoupled(K,T,L,0.01*np.eye(2,2))

In [6]:
KY,B,D = alg.Control_Decentral(K,T,L)
# 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_n9/Masterthesis_Models_mimo_0closedloop.fmu")
sim.setOperationMode('FMU for ModelExchange')
# Parameter Values
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 
        # System Gain
        params.update({"fmu.num["+str(outputs)+","+str(inputs)+",1]": k[outputs-1][inputs-1]})
        # System Lag
        for order in range(0,9):
            params.update({"fmu.den["+str(outputs)+","+str(inputs)+","+str(order+1)+"]": t[outputs-1][inputs-1][order]})
        # System Delay
        params.update({"fmu.delay["+str(outputs)+","+str(inputs)+"]": l[outputs-1][inputs-1]})
# Set Parameter Values
sim.set(params)
#sim.showParameterDialog()

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

p = bk.figure(title = "Test")
p.line(res["time"],res["fmu.y[1]"], color = "blue")
p.line(res["time"],res["fmu.y[2]"], color = "red")
bk.show(p)

In [7]:
# Load the closed loop model
sim.clear()
sim.loadModel("C:/Users/juliu/Documents/Thesis/Modelica/FMU/2_2_n9/Masterthesis_Models_mimo_0closedloop.fmu")
sim.setOperationMode('FMU for ModelExchange')
# Parameter Values
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 
        # System Gain
        params.update({"fmu.num["+str(outputs)+","+str(inputs)+",1]": k[outputs-1][inputs-1]})
        # System Lag
        for order in range(0,9):
            params.update({"fmu.den["+str(outputs)+","+str(inputs)+","+str(order+1)+"]": t[outputs-1][inputs-1][order]})
        # System Delay
        params.update({"fmu.delay["+str(outputs)+","+str(inputs)+"]": l[outputs-1][inputs-1]})
# Set Parameter Values
sim.set(params)
#sim.showParameterDialog()
sim.experimental_findSteadyState()
sim.analyser_getStateSpaceForm()


Out[7]:
{'A': array([[ -1.51036183e+00,  -9.90577679e-01,  -3.68619582e-01, ...,
           1.33332374e-06,   2.54109884e-16,   1.69406589e-16],
        [  1.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   1.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
        ..., 
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00]]),
 'B': array([[  2.94824260e-08,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   4.53318708e-08],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  2.28924256e-07,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  3.38813179e-16,   3.30273071e-07],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  0.00000000e+00,   0.00000000e+00],
        [  1.15007373e-03,   0.00000000e+00],
        [  0.00000000e+00,   1.00000000e-10],
        [  1.00000000e-10,   0.00000000e+00],
        [  0.00000000e+00,   1.82577921e-03]]),
 'C': array([[ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. , 10. ,  0. ,  0. ,  0. ,
          0. ,  0. ,  0. ,  0. ,  0.5,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,
          0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,
          0. ,  0. ,  0. ],
        [ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,
          0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,
          0. ,  5. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  8. ,  0. ,
          0. ,  0. ,  0. ]]),
 'D': array([[ 0.,  0.],
        [ 0.,  0.]]),
 'names_dxdt': ['fmu.System.miso[1].transferFunction[1].der(x_scaled[1])',
  'fmu.System.miso[1].transferFunction[1].der(x_scaled[2])',
  'fmu.System.miso[1].transferFunction[1].der(x_scaled[3])',
  'fmu.System.miso[1].transferFunction[1].der(x_scaled[4])',
  'fmu.System.miso[1].transferFunction[1].der(x_scaled[5])',
  'fmu.System.miso[1].transferFunction[1].der(x_scaled[6])',
  'fmu.System.miso[1].transferFunction[1].der(x_scaled[7])',
  'fmu.System.miso[1].transferFunction[1].der(x_scaled[8])',
  'fmu.System.miso[1].transferFunction[2].der(x_scaled[1])',
  'fmu.System.miso[1].transferFunction[2].der(x_scaled[2])',
  'fmu.System.miso[1].transferFunction[2].der(x_scaled[3])',
  'fmu.System.miso[1].transferFunction[2].der(x_scaled[4])',
  'fmu.System.miso[1].transferFunction[2].der(x_scaled[5])',
  'fmu.System.miso[1].transferFunction[2].der(x_scaled[6])',
  'fmu.System.miso[1].transferFunction[2].der(x_scaled[7])',
  'fmu.System.miso[1].transferFunction[2].der(x_scaled[8])',
  'fmu.System.miso[2].transferFunction[1].der(x_scaled[1])',
  'fmu.System.miso[2].transferFunction[1].der(x_scaled[2])',
  'fmu.System.miso[2].transferFunction[1].der(x_scaled[3])',
  'fmu.System.miso[2].transferFunction[1].der(x_scaled[4])',
  'fmu.System.miso[2].transferFunction[1].der(x_scaled[5])',
  'fmu.System.miso[2].transferFunction[1].der(x_scaled[6])',
  'fmu.System.miso[2].transferFunction[1].der(x_scaled[7])',
  'fmu.System.miso[2].transferFunction[1].der(x_scaled[8])',
  'fmu.System.miso[2].transferFunction[2].der(x_scaled[1])',
  'fmu.System.miso[2].transferFunction[2].der(x_scaled[2])',
  'fmu.System.miso[2].transferFunction[2].der(x_scaled[3])',
  'fmu.System.miso[2].transferFunction[2].der(x_scaled[4])',
  'fmu.System.miso[2].transferFunction[2].der(x_scaled[5])',
  'fmu.System.miso[2].transferFunction[2].der(x_scaled[6])',
  'fmu.System.miso[2].transferFunction[2].der(x_scaled[7])',
  'fmu.System.miso[2].transferFunction[2].der(x_scaled[8])',
  'fmu.Decentral_Controller.pi[1].pIController_sw[1].der(integral)',
  'fmu.Decentral_Controller.pi[1].pIController_sw[2].der(integral)',
  'fmu.Decentral_Controller.pi[2].pIController_sw[1].der(integral)',
  'fmu.Decentral_Controller.pi[2].pIController_sw[2].der(integral)'],
 'names_u': ['fmu.u[1]', 'fmu.u[2]'],
 'names_x': ['fmu.System.miso[1].transferFunction[1].x_scaled[1]',
  'fmu.System.miso[1].transferFunction[1].x_scaled[2]',
  'fmu.System.miso[1].transferFunction[1].x_scaled[3]',
  'fmu.System.miso[1].transferFunction[1].x_scaled[4]',
  'fmu.System.miso[1].transferFunction[1].x_scaled[5]',
  'fmu.System.miso[1].transferFunction[1].x_scaled[6]',
  'fmu.System.miso[1].transferFunction[1].x_scaled[7]',
  'fmu.System.miso[1].transferFunction[1].x_scaled[8]',
  'fmu.System.miso[1].transferFunction[2].x_scaled[1]',
  'fmu.System.miso[1].transferFunction[2].x_scaled[2]',
  'fmu.System.miso[1].transferFunction[2].x_scaled[3]',
  'fmu.System.miso[1].transferFunction[2].x_scaled[4]',
  'fmu.System.miso[1].transferFunction[2].x_scaled[5]',
  'fmu.System.miso[1].transferFunction[2].x_scaled[6]',
  'fmu.System.miso[1].transferFunction[2].x_scaled[7]',
  'fmu.System.miso[1].transferFunction[2].x_scaled[8]',
  'fmu.System.miso[2].transferFunction[1].x_scaled[1]',
  'fmu.System.miso[2].transferFunction[1].x_scaled[2]',
  'fmu.System.miso[2].transferFunction[1].x_scaled[3]',
  'fmu.System.miso[2].transferFunction[1].x_scaled[4]',
  'fmu.System.miso[2].transferFunction[1].x_scaled[5]',
  'fmu.System.miso[2].transferFunction[1].x_scaled[6]',
  'fmu.System.miso[2].transferFunction[1].x_scaled[7]',
  'fmu.System.miso[2].transferFunction[1].x_scaled[8]',
  'fmu.System.miso[2].transferFunction[2].x_scaled[1]',
  'fmu.System.miso[2].transferFunction[2].x_scaled[2]',
  'fmu.System.miso[2].transferFunction[2].x_scaled[3]',
  'fmu.System.miso[2].transferFunction[2].x_scaled[4]',
  'fmu.System.miso[2].transferFunction[2].x_scaled[5]',
  'fmu.System.miso[2].transferFunction[2].x_scaled[6]',
  'fmu.System.miso[2].transferFunction[2].x_scaled[7]',
  'fmu.System.miso[2].transferFunction[2].x_scaled[8]',
  'fmu.Decentral_Controller.pi[1].pIController_sw[1].integral',
  'fmu.Decentral_Controller.pi[1].pIController_sw[2].integral',
  'fmu.Decentral_Controller.pi[2].pIController_sw[1].integral',
  'fmu.Decentral_Controller.pi[2].pIController_sw[2].integral'],
 'names_y': ['fmu.y[1]', 'fmu.y[2]']}

In [24]:
# Load a Model
sim = mb.Simulator()
sim.clear()
sim.loadModel("C:/Users/juliu/Documents/Thesis/Modelica/FMU/2_2_n9/Masterthesis_Models_mimo_0processmodel.fmu")
sim.setOperationMode('FMU for ModelExchange')

# Show log window
sim.showLogWindow()

# 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
        for order in range(0,9):
            params.update({"fmu.den["+str(outputs)+","+str(inputs)+","+str(order+1)+"]": t[outputs-1][inputs-1][order]})
        #params.update({"fmu.den["+str(outputs)+","+str(inputs)+",9]": 1})
        # System Delay
        params.update({"fmu.delay["+str(outputs)+","+str(inputs)+"]": l[outputs-1][inputs-1]})

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

ss = sim.analyser_getStateSpaceForm()

In [27]:
A = ss['A']
B = ss['B']
C = ss['C']
D = ss['D']
I = np.eye(A.shape[0])

# Make Transfer function Matrix for frequency
G = np.dot(np.dot(C,np.linalg.inv(0*1j*I-A)),B)+D
S = np.inv()
G = np.absolute(G)
U,V,W = np.linalg.svd(G)

In [29]:



Out[29]:
array([[-0.74283446, -0.66947514],
       [-0.66947514,  0.74283446]])

In [ ]: