* cost_function edits, const modeltype, #55

pyDMPC/ControlFramework/Init.py

@@ -1,9 +1,9 @@
 # Global paths
-glob_lib_paths = [r'C:\Git\pyDMPC\pyDMPC\ModelicaModels\ModelicaModels',
-             r'C:\Git\modelica-buildings\Buildings',
-             r'C:\Git\AixLib\AixLib']
-glob_res_path = r'C:\TEMP\Dymola'
-glob_dym_path = r'C:\Program Files\Dymola 2018 FD01\Modelica\Library\python_interface\dymola.egg'
+glob_lib_paths = [r'C:\mst\pyDMPC\pyDMPC\ModelicaModels\ModelicaModels',
+             r'C:\mst\modelica-buildings\Buildings',
+             r'C:\mst\AixLib\AixLib']
+glob_res_path = r'C:\mst\dymola'
+glob_dym_path = r'C:\Program Files (x86)\Dymola 2018\Modelica\Library\python_interface\dymola.egg'
 
 # Working directory
 import time
@@ -48,6 +48,7 @@
 
 # Modifiers
 cost_fac = []
+factors = []
 
 # Variation
 min_var = []
@@ -63,57 +64,61 @@
 
 # Subsystems
 sys_id.append(0)
-name.append("Field")
-model_type.append("Modelica")
+name.append("Building")
+model_type.append("Linear")
 ups_neigh.append(None)
 downs_neigh.append(1)
-input_names.append(["returnTemperature.T"])
-input_variables.append(["external"])
-inputs.append([])
-output_names.append(["returnTemperature.T"])
+input_names.append(["supplyTemperature.T"])
+input_variables.append([r"variation.table[1,2]"])
+inputs.append(range(280,290,5))
+output_names.append(["returnTemperature"])
 set_points.append([287])
-state_var_names.append(["supplyTemperature.T"])
-model_state_var_names.append(["vol.T_start"])
+state_var_names.append(["const.y"])
+model_state_var_names.append(["const.k"])   #"const.k"
 start.append(0.)
-stop.append(3600.0*24*365.25*3)
-incr.append(3600.)
-opt_time.append(0)
+stop.append(7200.)
+incr.append(10.)
+opt_time.append(3600)
 samp_time.append(10)
 lib_paths.append(glob_lib_paths)
 res_path.append(glob_res_path + "\\" + name_wkdir)
 dym_path.append(glob_dym_path)
-mod_path.append(r'ModelicaModels.SubsystemModels.DetailedModels.Geo.Field')
+mod_path.append(r'ModelicaModels.SubsystemModels.DetailedModels.Geo.Building')
 command_names.append(["heatShare"])
 command_variables.append(["decisionVariables.table[1,2]"])
 commands.append(range(0,105,5))
-traj_points.append(range(278,310,1))
-traj_var.append(["supplyTemperature.T"])
-cost_fac.append([0.0, 0.0, 1.0, 0.1])
+traj_points.append([])
+traj_var.append([])
+cost_fac.append([-1.0, 0.0, 1000.0, -1000.0, 0.0, 0.0, 0.0, 0.0])
+factors.append([1, -60./4.18/8./100.])
 
 sys_id.append(1)
-name.append("Building")
-model_type.append("Modelica")
+name.append("Field")
+model_type.append("const")
 ups_neigh.append(0)
 downs_neigh.append(None)
-input_names.append(["supplyTemperature.T"])
-input_variables.append([r"variation.table[1,2]"])
-inputs.append(range(280,310,10))
-output_names.append(["returnTemperature"])
-set_points.append([287])
-state_var_names.append(["sine.y"])
-model_state_var_names.append(["const.k"])
+input_names.append(["returnTemperature.T"])
+input_variables.append(["external"])
+inputs.append([])
+output_names.append(["movMea.y"])
+set_points.append([285.65])
+state_var_names.append(["supplyTemperature.T"])
+model_state_var_names.append(["vol.T_start"])
 start.append(0.)
-stop.append(7200.)
-incr.append(10.)
-opt_time.append(0)
+stop.append(3600.0*24*365.25*3)
+incr.append(3600.)
+opt_time.append(86400)
 samp_time.append(10)
 lib_paths.append(glob_lib_paths)
 res_path.append(glob_res_path + "\\" + name_wkdir)
 dym_path.append(glob_dym_path)
-mod_path.append(r'ModelicaModels.SubsystemModels.DetailedModels.Geo.Building')
-command_names.append([])
+mod_path.append(r'ModelicaModels.SubsystemModels.DetailedModels.Geo.Field')
+command_names.append(["traj"])
 command_variables.append(["decisionVariables.table[1,2]"])
 commands.append(range(0,105,5))
-traj_points.append([])
-traj_var.append([])
-cost_fac.append([-0.01, 1.0, 0.0, 0.1])
+traj_points.append(range(280,290,1))
+traj_var.append(["supplyTemperature.T"])
+cost_fac.append([-1.0, 0.0, 10.0, -10.0, 0.0, 0.0, 0.0, 0.0 ])   
+factors.append([0,0])
+#Legend: cost_fac.append([real_cost, cost_downstr_neigh, cost_dev_setpoint_penalty, cost_dev_setpoint_reward, cost_integ_penalty, cost_integ_reward, cost_diff_penalty, cost_diff_reward])
+

pyDMPC/ControlFramework/Modeling.py

@@ -253,6 +253,36 @@ def get_results(self, name):
     def predict(self):
         self.simulate()
         self.get_outputs()
+
+class ConstMod(Model):
+
+    dymola = None
+    lib_paths = Init.glob_lib_paths
+    dym_path = Init.glob_dym_path
+
+    def __init__(self, sys_id):
+        super().__init__(sys_id)
+
+    def get_outputs(self):
+        self.states.outputs = []
+
+        if self.states.output_names is not None:
+            for nam in self.states.output_names:
+                self.states.outputs.append(self.get_results(nam))
+
+    def get_results(self, name):
+        from modelicares import SimRes
+        import os
+        # Get the simulation result
+        sim = SimRes(os.path.join(self.paths.res_path, 'dsres.mat'))
+
+        arr = sim[name].values()
+        results = arr.tolist()
+
+        return results
+
+    def predict(self):
+        self.get_outputs()
 
 class SciMod(Model):
     def __init__(self, sys_id):

pyDMPC/ControlFramework/Subsystem.py

@@ -2,6 +2,7 @@
 import Modeling
 import System
 import Time
+import numpy as np
 
 class Subsystem:
 
@@ -97,6 +98,8 @@ def __init__(self, sys_id):
         self.traj_var = Init.traj_var[sys_id]
         self.traj_points = Init.traj_points[sys_id]
         self.counter = 0
+        self.setpoint_prev = 0
+        self.phase = 0
 
     def prepare_model(self):
         """Prepares the model of a subsystem according to the subsystem's model
@@ -119,6 +122,8 @@ def prepare_model(self):
             model.load_mod()
         elif self.model_type == "Linear":
             model = Modeling.LinMod(self.sys_id)
+        elif self.model_type == "const":
+            model = Modeling.ConstMod(self.sys_id)
         else:
             model = Modeling.FuzMod(self.sys_id)
         return model
@@ -200,6 +205,10 @@ def interp_minimize(self, interp):
         if self.traj_var != []:
             traj_costs = []
             traj = self.model.get_results(self.traj_var[0])
+            if self.setpoint_send:
+                self.setpoint_prev = self.setpoint_send
+            else:
+                self.setpoint_prev = 285.65
             self.setpoint_send = traj[10]
 
         else:
@@ -233,7 +242,6 @@ def interp_minimize(self, interp):
             else:
                 self.coup_vars_send = opt_outputs[0]
 
-
     def calc_cost(self, command, outputs):
         import scipy.interpolate
 
@@ -257,31 +265,50 @@ def calc_cost(self, command, outputs):
             setpoint = self.setpoint_rec
         else:
             setpoint = self.model.states.set_points[0]
-
-        if self.model.states.set_points != []:
-            self.err_prop += (outputs -setpoint)
+
+        if self.setpoint_send:
+            self.phase = self.setpoint_prev - self.setpoint_send    #if > 0: field for heating the building, if < 0: field for cooling the building 
+
+        #Regelung in Anlehnung an PID Regelverhalten
+        #Proportionalanteil
+        self.err_prop = outputs - setpoint                          #Deviation from setpoint (proportional & integral)
+
+        if self.phase > 0:#if field cools off
+            if self.err_prop < 0:                                  
+                cost += self.cost_fac[2]*self.err_prop              #Penalty Deviation (proportional)
+            else:
+                cost += self.cost_fac[3]*(setpoint - outputs)           #Reward Deviation (proportional)
+        else:#if field heats up
             if self.err_prop > 0:
-                cost += (self.cost_fac[2] * (self.err_prop)**2)
+                cost += self.cost_fac[2]*self.err_prop              #Penalty Deviation (proportional)
             else:
-                cost += (self.cost_fac[3] * (self.err_prop)**2)
-
+                cost += self.cost_fac[3]*(setpoint-outputs)           #Reward Deviation (proportional)
+
+        #Integralteil
+        self.err_integ += outputs - setpoint
+
         if self.cost_rec != []:    
-            self.err_integ += (outputs - setpoint)
-            if self.err_integ > 0:
-                cost += self.cost_fac[3] * self.err_integ           #Integral penalization
+            if self.phase > 0:                                      #if field cools off
+                if self.err_integ > 0:
+                    cost += self.cost_fac[4] * self.err_integ       #Integral reward
+                else:
+                    cost += self.cost_fac[3] * (-(self.err_integ))       #Integral penalty
             else:
-                cost += self.cost_fac[4] * self.err_integ           #Integral reward
-
-            self.counter += 1
-            self.err_prev = self.err_curr
-            self.err_curr = outputs - setpoint
-
-            if self.counter > 1:
-                self.err_diff = self.err_curr - self.err_prev
-                if self.err_diff > 0:
-                    cost += self.cost_fac[5] * self.err_diff        #Differential penalization
+                if self.err_integ > 0:
+                    cost += self.cost_fac[3] * self.err_integ       #Integral penalization
                 else:
-                    cost += self.cost_fac[6] * self.err_diff        #Differential reward
+                    cost += self.cost_fac[4] * (-(self.err_integ))       #Integral reward
+
+#        self.counter += 1
+#        self.err_prev = self.err_curr
+#        self.err_curr = outputs - setpoint
+#            
+#        if self.counter > 1:
+#            self.err_diff = self.err_curr - self.err_prev
+#            if self.err_diff > 0:
+#                cost += self.cost_fac[5] * self.err_diff        #Differential penalization
+#            else:
+#                cost += self.cost_fac[6] * self.err_diff        #Differential reward
 
         return cost
 

pyDMPC/ControlFramework/System.py

@@ -62,13 +62,17 @@ def gen_subsys(self):
             The created subsystem objects
         """
         import os
+        import shutil
         subsystems = []
 
         for i in range(self.amo_subsys):
             subsystems.append(Subsystem.Subsystem(i))
 
         for sys in subsystems:
             os.mkdir(self.wkdir + "\\" + sys.name)
+            if sys.model_type == "const":
+                shutil.copyfile(Init.glob_res_path + "\\dsres.mat", 
+                        sys.model.paths.res_path + "\\dsres.mat")
 
         return subsystems
 

pyDMPC/ModelicaModels/ModelicaModels/ControlledSystems/GeothermalField.mo

@@ -2,20 +2,19 @@ within ModelicaModels.ControlledSystems;
 model GeothermalField
 
   extends ModelicaModels.Subsystems.Geo.BaseClasses.FieldBaseClass(
-    vol(T_start=287.15),
-    pump(T_start=287.15),
-    vol1(V=2, T_start=287.15),
+    vol(T_start=282.15),
+    pump(T_start=282.15),
+    vol1(V=2, T_start=282.15),
     fixedTemperature(T=285.65),
-    pressurePoint(use_T_in=false, T=285.65),
+    pressurePoint(use_T_in=false, T=282.15),
     supplyTemperature(T(fixed=true, start=285.65)),
-    returnTemperature(T(fixed=true, start=285.65)));
+    returnTemperature(T(fixed=true, start=285.65)),
+    thermalConductor(G=1000));
   Modelica.Thermal.HeatTransfer.Sources.PrescribedHeatFlow prescribedHeatFlow
     annotation (Placement(transformation(extent={{-12,-12},{12,12}},
         rotation=0,
         origin={-6,42})));
-  Modelica.Blocks.Sources.Constant
-                               const(k=-60000)
-    annotation (Placement(transformation(extent={{-100,38},{-80,60}})));
+  Modelica.Blocks.Sources.Constant Error;
   Modelica.Blocks.Math.Product product
     annotation (Placement(transformation(extent={{-62,32},{-42,52}})));
   Modelica.Blocks.Interfaces.RealInput heatShare
@@ -30,11 +29,11 @@ model GeothermalField
         origin={-74,20})));
   Modelica.Thermal.HeatTransfer.Celsius.FromKelvin fromKelvin
     annotation (Placement(transformation(extent={{-74,-60},{-54,-40}})));
+  Modelica.Blocks.Sources.Constant const(k=-60000)
+    annotation (Placement(transformation(extent={{-100,38},{-80,58}})));
 equation
   connect(prescribedHeatFlow.port, vol1.heatPort)
     annotation (Line(points={{6,42},{20,42}}, color={191,0,0}));
-  connect(const.y, product.u1) annotation (Line(points={{-79,49},{-74,49},{-74,48},
-          {-64,48}}, color={0,0,127}));
   connect(heatShare, percent.u)
     annotation (Line(points={{-100,2},{-74,2},{-74,12.8}}, color={0,0,127}));
   connect(percent.y, product.u2)
@@ -43,6 +42,8 @@ equation
     annotation (Line(points={{-41,42},{-18,42}}, color={0,0,127}));
   connect(traj, fromKelvin.Kelvin)
     annotation (Line(points={{-100,-50},{-76,-50}}, color={0,0,127}));
+  connect(const.y, product.u1)
+    annotation (Line(points={{-79,48},{-64,48}}, color={0,0,127}));
   annotation (Icon(coordinateSystem(preserveAspectRatio=false)), Diagram(
         coordinateSystem(preserveAspectRatio=false)));
 end GeothermalField;

pyDMPC/ModelicaModels/ModelicaModels/DataBase.mo

@@ -4,7 +4,7 @@ package DataBase "Contians records for the various case studies"
 
     record GeoRecord
 
-      parameter Modelica.SIunits.MassFlowRate m_flow_tot = 16.0
+      parameter Modelica.SIunits.MassFlowRate m_flow_tot = 8.0
        "The total mass flow rate circulating through field and
    building";
 

pyDMPC/ModelicaModels/ModelicaModels/SubsystemModels/DetailedModels/Geo/Building.mo

@@ -4,7 +4,9 @@ model Building "Simplified building model"
   extends ModelicaModels.SubsystemModels.DetailedModels.Geo.GeoCommunicationBaseClass(variation(
         table=[0,293]));
 
-  extends ModelicaModels.Subsystems.Geo.BaseClasses.BuildingBaseClass(const(k=-10000));
+  extends ModelicaModels.Subsystems.Geo.BaseClasses.BuildingBaseClass(const(k=-10000),
+    vol1(T_start=281.15),
+    bou(T=281.15));
 
 equation
   connect(variation.y[1], boundary.T_in) annotation (Line(points={{-79,90},{-68,

pyDMPC/ModelicaModels/ModelicaModels/SubsystemModels/DetailedModels/Geo/Field.mo

@@ -3,23 +3,23 @@ model Field "Simplified model of geothermal field"
 
   extends
     ModelicaModels.SubsystemModels.DetailedModels.Geo.GeoCommunicationBaseClass(variation(
-        extrapolation=Modelica.Blocks.Types.Extrapolation.Periodic, table=[0,10000;
-          2635200,12000; 5270400,9000; 7905600,3000; 10540800,-5000; 13176000,-12000;
-          15811200,-14000; 18446400,-13000; 21081600,-6000; 23716800,1000; 26352000,
-          4000; 28987200,8000]),                          decisionVariables(
+        extrapolation=Modelica.Blocks.Types.Extrapolation.Periodic, table=[0,
+          10000; 2635200,12000; 5270400,9000; 7905600,3000; 10540800,-5000;
+          13176000,-12000; 15811200,-14000; 18446400,-13000; 21081600,-6000;
+          23716800,1000; 26352000,4000; 28987200,8000]),  decisionVariables(
         extrapolation=Modelica.Blocks.Types.Extrapolation.Periodic, table=[0.0,30]),
     percent(k=0.01));
 
   extends ModelicaModels.Subsystems.Geo.BaseClasses.FieldBaseClass(
       returnTemperature(T(start=285.65, displayUnit="K")),
-      vol1(V=2, T_start=287.15),
+      vol1(V=2, T_start=281.15),
       supplyTemperature(T(start=285.65, displayUnit="K")),
-      pressurePoint(T(start=285.65, fixed=true)),
-      pump(addPowerToMedium=false, T_start=287.15),
+      pressurePoint(T(start=281.15, fixed=true)),
+      pump(addPowerToMedium=false, T_start=281.15),
     movMea(delta=94672800),
-    vol(V=6000, T_start=287.15),
+    vol(V=6000, T_start=281.15),
     fixedTemperature(T=285.65),
-    thermalConductor(G=50));
+    thermalConductor(G=1000));
 
   Modelica.Blocks.Math.Product product1
     annotation (Placement(transformation(extent={{-78,0},{-58,20}})));