changed the progressbar implementation and updated unittests to run t…

…he examples

doc/examples/quickstart.py

@@ -0,0 +1,189 @@
+#!/usr/bin/env python
+################################################################################
+#    Copyright 2015 Brecht Baeten
+#    This file is part of mpcpy.
+#
+#    mpcpy is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation, either version 3 of the License, or
+#    (at your option) any later version.
+#
+#    mpcpy is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    You should have received a copy of the GNU General Public License
+#    along with mpcpy.  If not, see <http://www.gnu.org/licenses/>.
+################################################################################
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import mpcpy
+import pyomo.environ as pyomo
+
+
+# Define an emulator class
+class Emulator(mpcpy.Emulator):
+    """
+    A custom system emulator
+    """
+    def simulate(self, starttime, stoptime, input):
+        dt = 1
+        time = np.arange(starttime, stoptime+dt, dt, dtype=np.float)
+
+        # initialize
+        x = np.ones_like(time)*self.res['x'][-1]
+
+        # interpolate inputs
+        u = np.interp(time, input['time'], input['u'])
+        d = np.interp(time, input['time'], input['d'])
+
+        # perform simulation
+        for i, t in enumerate(time[:-1]):
+            # dx/dt = A*x + d + u
+            x[i+1] = x[i] + (self.parameters['A']*x[i] + d[i] + u[i])*dt
+
+        # create and return a results dict
+        res = {
+            'time': time,
+            'x': x,
+            'd': d,
+            'u': u,
+        }
+
+        return res
+
+
+# Define a control class
+class SetpointControl(mpcpy.Control):
+    """
+    A control to keep the state as close to a set point as possible
+    """
+    def formulation(self):
+        # create a pyomo model
+        model = pyomo.AbstractModel()
+
+        model.i = pyomo.Set()
+        model.ip = pyomo.Set()
+
+        model.time = pyomo.Param(model.ip)
+        model.d = pyomo.Param(model.ip, initialize=0.)
+        model.x = pyomo.Var(model.ip, domain=pyomo.Reals, initialize=0.)
+        model.u = pyomo.Var(model.ip, domain=pyomo.NonNegativeReals, bounds=(0., 1.), initialize=0.)
+
+        model.x0 = pyomo.Param(initialize=0.)
+
+        model.initialcondition = pyomo.Constraint(
+            rule=lambda model: model.x[0] == model.x0
+        )
+
+        model.constraint = pyomo.Constraint(
+            model.i,
+            rule=lambda model, i: (model.x[i+1]-model.x[i])/(model.time[i+1]-model.time[i]) ==
+                          self.parameters['A']*model.x[i] + model.d[i] + model.u[i]
+        )
+
+        model.objective = pyomo.Objective(
+            rule=lambda model: sum((model.x[i]-self.parameters['set'])**2 for i in model.i)
+        )
+
+        # store the model inside the object
+        self.model = model
+
+    def solution(self, sta, pre):
+        # create data and instantiate the pyomo model
+        ip = np.arange(len(pre['time']))
+        data = {
+            None: {
+                'i': {None: ip[:-1]},
+                'ip': {None: ip},
+                'time': {(i,): v for i, v in enumerate(pre['time'])},
+                'x0': {None: sta['x']},
+                'd': {(i,): pre['d'][i] for i in ip},
+            }
+        }
+
+        instance = self.model.create_instance(data)
+
+        # solve and return the control inputs
+        optimizer = pyomo.SolverFactory('ipopt')
+        results = optimizer.solve(instance)
+
+        sol = {
+            'time': np.array([pyomo.value(instance.time[i]) for i in instance.ip]),
+            'x': np.array([pyomo.value(instance.x[i]) for i in instance.ip]),
+            'u': np.array([pyomo.value(instance.u[i]) for i in instance.ip]),
+            'd': np.array([pyomo.value(instance.d[i]) for i in instance.ip]),
+        }
+
+        return sol
+
+
+# Define a state estimation class
+class StateestimationPerfect(mpcpy.Stateestimation):
+    """
+    Perfect state estimation
+    """
+    def stateestimation(self, time):
+        return {'x': np.interp(time, self.emulator.res['time'], self.emulator.res['x'])}
+
+
+# instantiate the emulator
+emulator = Emulator(['u', 'd'], parameters={'A': -0.2}, initial_conditions={'x': 0})
+
+# test the emulator with some random data
+time = np.arange(0., 1001., 10.)
+np.random.seed(0)
+d = np.random.random(len(time)) - 0.5
+u = 1.0*np.ones(len(time))
+
+emulator.initialize()
+res = emulator(time, {'time': time, 'd': d, 'u': u})
+print(res)
+
+
+# create a disturbances object
+time = np.arange(0., 1001., 10.)
+d = 0.5*np.sin(2*np.pi*time/1000)
+disturbances = mpcpy.Disturbances({'time': time, 'd': d})
+
+bcs = disturbances(np.array([0, 20, 40, 60, 100]))
+print(bcs)
+
+
+# create a stateestimation object
+stateestimation = StateestimationPerfect(emulator)
+sta = stateestimation(0)     
+print(sta)        
+
+
+# create a prediction object
+prediction = mpcpy.Prediction(disturbances)
+pre = prediction(np.array([0, 20, 40, 60, 100]))
+print(pre)   
+
+
+# create a control object and mpc object
+control = SetpointControl(stateestimation, prediction, parameters={'A': -0.2, 'set': 3.0},
+                          horizon=100., timestep=10., receding=10.)
+mpc = mpcpy.MPC(emulator, control, disturbances, emulationtime=1000, resulttimestep=10)
+
+# run the mpc
+res = mpc(verbose=1)
+print(res)
+
+
+# plot results
+fig, ax = plt.subplots(2, 1)
+ax[0].plot(res['time'], res['u'])
+ax[0].set_ylabel('u')
+
+ax[1].plot(res['time'], res['x'])
+ax[1].set_xlabel('time')
+ax[1].set_ylabel('x')
+
+
+if __name__ == '__main__':
+    plt.show()

doc/source/mpc.rst

@@ -3,4 +3,4 @@ MPC
 
 .. autoclass:: mpcpy.MPC
    :members: 
-
+   :special-members: __call__

examples/quickstart.py

@@ -171,20 +171,19 @@ def stateestimation(self, time):
 mpc = mpcpy.MPC(emulator, control, disturbances, emulationtime=1000, resulttimestep=10)
 
 # run the mpc
-res = mpc()
+res = mpc(verbose=1)
 print(res)
 
 
-def plot_result():
-    fig, ax = plt.subplots(2, 1)
-    ax[0].plot(res['time'], res['u'])
-    ax[0].set_ylabel('u')
+# plot results
+fig, ax = plt.subplots(2, 1)
+ax[0].plot(res['time'], res['u'])
+ax[0].set_ylabel('u')
 
-    ax[1].plot(res['time'], res['x'])
-    ax[1].set_xlabel('time')
-    ax[1].set_ylabel('x')
+ax[1].plot(res['time'], res['x'])
+ax[1].set_xlabel('time')
+ax[1].set_ylabel('x')
 
 
 if __name__ == '__main__':
-    plot_result()
     plt.show()

examples/simple_space_heating_mpc.py

@@ -262,7 +262,7 @@ def solution(self, sta, pre):
 
 # MPC
 mpc = mpcpy.MPC(emulator, control, disturbances, emulationtime=1*24*3600., resulttimestep=60)
-res = mpc()
+res = mpc(verbose=1)
 
 
 # Plot results
@@ -289,7 +289,7 @@ def solution(self, sta, pre):
 def_control = LinearProgram(def_stateestimation, prediction,
                             parameters=control_parameters, horizon=24*3600., timestep=3600.)
 def_mpc = mpcpy.MPC(def_emulator, def_control, disturbances, emulationtime=1*24*3600., resulttimestep=60)
-def_res = def_mpc()
+def_res = def_mpc(verbose=1)
 
 fix, ax = plt.subplots(2, 1)
 ax[0].plot(def_res['time']/3600, def_res['Q_flow_hp'], 'k', label='hp')
@@ -308,4 +308,5 @@ def solution(self, sta, pre):
 ax[1].legend(loc='lower right')
 
 
-plt.show()
+if __name__ == '__main__':
+    plt.show()

mpcpy/mpc.py

@@ -74,10 +74,15 @@ def nextstepcalculator(controlsolution):
         self.res = {}
         self.appendres = {}
 
-    def __call__(self):
+    def __call__(self, verbose=0):
         """
         Runs the mpc simulation
-         
+
+        Parameters
+        ----------
+        verbose: optional, int
+            Controls the amount of print output
+
         Returns
         -------
         dict
@@ -88,19 +93,17 @@ def __call__(self):
         # initialize the emulator
         self.emulator.initialize()
         starttime = 0
-
-
+
         if self.plotfunction:
             (fig,ax,pl) = self.plotfunction()
 
         # prepare a progress bar
-        barwidth = 80
+        barwidth = 80-2
         barvalue = 0
-        print('Run MPC %s |' %(' '*(barwidth-10)))
-        # sys.stdout.write('[%s]\n' % (' ' * barwidth))
-        # sys.stdout.flush()
-        # sys.stdout.write('\b' * (barwidth+1))
-
+        if verbose > 0:
+            print('Running MPC')
+            print('[' + (' '*barwidth) + ']', end='')
+
         while starttime < self.emulationtime:
 
             # calculate control signals for the control horizon
@@ -147,21 +150,18 @@ def __call__(self):
             starttime = self.emulator.res['time'][-1]
 
             # update the progress bar
-            if starttime/self.emulationtime*barwidth >= barvalue:
-                addbars = int(round(starttime/self.emulationtime*barwidth-barvalue))
-                sys.stdout.write(addbars*'-')
-                sys.stdout.flush()
-                barvalue += addbars
-
+            if verbose > 0:
+                if starttime/self.emulationtime*barwidth >= barvalue:
+                    barvalue += int(round(starttime/self.emulationtime*barwidth-barvalue))
+                    print('\r[' + ('='*barvalue) + (' '*(barwidth-barvalue)) + ']', end='')
+
         # copy the results to a local res dictionary
         self.res.update(self.emulator.res)
 
         # interpolate the boundary conditions and add them to self.res
         self.res.update(self.disturbances(self.res['time']))
-
-        sys.stdout.write('  done')
-        sys.stdout.write("\n")
-        sys.stdout.flush()
+        if verbose > 0:
+            print(' done')
 
         return self.res
 

tests/examples.py

@@ -17,40 +17,18 @@
 #    along with mpcpy.  If not, see <http://www.gnu.org/licenses/>.
 ################################################################################
 
-import unittest
-import mpcpy
-import numpy as np
-import sys
-import os
-import subprocess
 
-# module path
-modulepath = os.path.abspath(os.path.dirname(sys.modules['mpcpy'].__file__))
-examplespath =  os.path.abspath(os.path.join(modulepath,'..','examples'))
+import unittest
 
-# define null file
-fnull = open(os.devnull, 'w')
 
 class TestExamples(unittest.TestCase):
 
     def test_simple_space_heating_mpc(self):
-
-        currentdir = os.getcwd()
-        os.chdir(examplespath)
-        filename = 'simple_space_heating_mpc'
-        p = subprocess.Popen(['jupyter', 'nbconvert', '--execute', '--to', 'notebook', '{}.ipynb'.format(filename)], stdout=fnull, stderr=subprocess.PIPE)
-        output, error = p.communicate()
-        os.remove('{}.nbconvert.ipynb'.format(filename))
-        os.chdir(currentdir)
-
-        self.assertEqual(p.returncode, 0, error)
+        from examples import simple_space_heating_mpc
 
     def test_quickstart(self):
-
         from examples import quickstart
 
-
-
 
 if __name__ == '__main__':
     unittest.main()