source: PyFMI/trunk/src/pyfmi/master.pyx

Last change on this file was 9970, checked in by Christian Andersson, 5 months ago

Fixed problem with step outside simulation region (ticket:5397) and added option for max step-size (ticket:5396)

File size: 76.4 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3
4# Copyright (C) 2010 Modelon AB
5#
6# This program is free software: you can redistribute it and/or modify
7# it under the terms of the GNU Lesser General Public License as published by
8# the Free Software Foundation, version 3 of the License.
9#
10# This program is distributed in the hope that it will be useful,
11# but WITHOUT ANY WARRANTY; without even the implied warranty of
12# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13# GNU Lesser General Public License for more details.
14#
15# You should have received a copy of the GNU Lesser General Public License
16# along with this program. If not, see <http://www.gnu.org/licenses/>.
17
18import pyfmi.fmi as fmi
19from pyfmi.fmi_algorithm_drivers import FMICSAlgOptions, AssimuloSimResult
20from pyfmi.common.algorithm_drivers import OptionBase, InvalidAlgorithmOptionException
21from pyfmi.common.io import ResultDymolaTextual, ResultHandlerFile, ResultHandlerDummy
22from pyfmi.common.core import TrajectoryLinearInterpolation
23from pyfmi.common.core import TrajectoryUserFunction
24
25from timeit import default_timer as timer
26
27import fnmatch
28import sys
29import re
30from collections import OrderedDict
31import time
32import numpy as np
33import warnings
34cimport numpy as np
35import scipy.sparse as sp
36import scipy.linalg as lin
37import scipy.sparse.linalg as splin
38import scipy.optimize as sopt
39import scipy.version
40
41from fmi cimport FMUModelCS2
42from cpython cimport bool
43cimport fmil_import as FMIL
44from fmi_util import Graph
45from cython.parallel import prange, parallel
46cimport openmp
47
48DEF SERIAL   = 0
49DEF PARALLEL = 1
50
51try:
52    from numpy.lib import NumpyVersion
53    USE_ROOT = NumpyVersion(scipy.version.version) >= "0.11.0"
54except ImportError: #Numpy version is < 1.9.0 so assume scipy version is the same
55    USE_ROOT = False
56
57cdef reset_models(list models):
58    for model in models:
59        model.reset()
60
61cdef perform_do_step(list models, dict time_spent, FMIL.fmi2_import_t** model_addresses, double cur_time, double step_size, bool new_step, int setting):
62    if setting == SERIAL:
63        perform_do_step_serial(models, time_spent, cur_time,  step_size,  new_step)
64    else:
65        perform_do_step_parallel(models, model_addresses, len(models), cur_time,  step_size,  new_step)
66
67cdef perform_do_step_serial(list models, dict time_spent, double cur_time, double step_size, bool new_step):
68    """
69    Perform a do step on all the models.
70    """
71    cdef double time_start = 0.0
72    cdef int status = 0
73   
74    for model in models: 
75        time_start = timer()
76        status = model.do_step(cur_time, step_size, new_step)
77        time_spent[model] += timer() - time_start
78       
79        if status != 0:
80            raise fmi.FMUException("The step failed for model %s at time %f. See the log for more information. Return flag %d."%(model.get_name(), cur_time, status))
81
82cdef perform_do_step_parallel(list models, FMIL.fmi2_import_t** model_addresses, int n, double cur_time, double step_size, int new_step):
83    """
84    Perform a do step on all the models.
85    """
86    cdef int i, status = 0
87    cdef int num_threads, id
88    cdef double time
89   
90    for i in prange(n, nogil=True, schedule="dynamic", chunksize=1): 
91    #for i in prange(n, nogil=True, schedule="dynamic"):
92        #num_threads = openmp.omp_get_num_threads()
93        #id = openmp.omp_get_thread_num()
94        #time = openmp.omp_get_wtime()
95       
96        status |= FMIL.fmi2_import_do_step(model_addresses[i], cur_time, step_size, new_step)
97        #time = openmp.omp_get_wtime() -time
98        #printf("Time: %f (%d), %d, %d, Elapsed Time: %f \n",cur_time, i, num_threads, id,  time)
99
100    if status != 0:
101        raise fmi.FMUException("The simulation failed. See the log for more information. Return flag %d."%status)
102   
103    #Update local times in models
104    for model in models:
105        model.time = cur_time + step_size
106
107cdef enter_initialization_mode(list models, double start_time, double final_time, object opts, dict time_spent):
108    cdef int status
109    for model in models:
110        time_start = timer()
111        model.setup_experiment(tolerance=opts["local_rtol"], start_time=start_time, stop_time=final_time)
112       
113        try:
114            status = model.enter_initialization_mode()
115            time_spent[model] += timer() - time_start
116        except fmi.FMUException:
117            print("The model, '" + model.get_name() + "' failed to enter initialization mode. ")
118       
119cdef exit_initialization_mode(list models, dict time_spent):
120    for model in models:
121        time_start = timer()
122        model.exit_initialization_mode()
123        time_spent[model] += timer() - time_start
124
125"""
126cdef perform_initialize(list models, double start_time, double final_time, object opts):
127    #
128    Initialize all the models.
129    #
130    for model in models:
131        model.setup_experiment(tolerance=opts["local_rtol"], start_time=start_time, stop_time=final_time)
132        model.initialize()
133"""
134
135#cdef get_fmu_states(list models):
136cdef get_fmu_states(list models, dict states_dict = None):
137    """
138    Get the FMU states for all the models
139    """
140    if states_dict is None:
141        return {model: model.get_fmu_state() for model in models}
142    else:
143        for model in states_dict.keys():
144            states_dict[model] = model.get_fmu_state(states_dict[model])
145        return states_dict
146   
147cdef set_fmu_states(dict states_dict):
148    """
149    Sets the FMU states for all the models
150    """
151    for model in states_dict.keys():
152        model.set_fmu_state(states_dict[model])
153       
154cdef free_fmu_states(dict states_dict):
155    """
156    Free the FMU states for all the models
157    """
158    for model in states_dict.keys():
159        model.free_fmu_state(states_dict[model])
160       
161cdef store_communication_point(object models_dict):
162    for model in models_dict.keys():
163        models_dict[model]["result"].integration_point()
164
165cdef finalize_result_objects(object models_dict):
166    for model in models_dict.keys():
167        models_dict[model]["result"].simulation_end()
168       
169def init_f(y, master):
170    y = y.reshape(-1, 1)
171    u = master.L.dot(y)
172    master.set_connection_inputs(u)
173    temp_y = y - master.get_connection_outputs()
174
175    return temp_y.flatten()
176   
177def init_f_block(ylocal, master, block):
178    y = np.zeros((master._len_outputs))
179    y[block["global_outputs_mask"]] = ylocal
180    y = y.reshape(-1,1)
181   
182    ytmp = np.zeros((master._len_outputs))
183    u = master.L.dot(y.reshape(-1,1))
184   
185    for model in block["inputs"].keys(): #Set the inputs
186        master.set_specific_connection_inputs(model, block["inputs_mask"][model], u)
187    for model in block["outputs"].keys(): #Get the outputs
188        master.get_specific_connection_outputs(model, block["outputs_mask"][model], ytmp)
189   
190    res = y - ytmp.reshape(-1,1)
191
192    return res.flatten()[block["global_outputs_mask"]]
193   
194def init_jac(y, master):
195    y = y.reshape(-1, 1)
196    u = master.L.dot(y)
197    master.set_connection_inputs(u)
198    D = master.compute_global_D()
199    DL = D.dot(master.L)
200    return np.eye(*DL.shape) - DL
201
202"""
203cdef np.ndarray matvec(object A, np.ndarray[np.float64_t, ndim=1, mode="c"] x):
204    cdef np.ndarray[dtype=np.int32_t, mode="c"] indices = A.indices
205    cdef np.ndarray[dtype=np.int32_t, mode="c"] indptr  = A.indptr
206    cdef np.ndarray[dtype=np.float64_t, mode="c"] data  = A.data
207    cdef np.ndarray b     = np.empty((x.size,1), dtype=np.float64)
208    cdef int i, j, ind
209    cdef double tmp
210   
211    for i in range(x.size):
212        #b[i] = np.sum(data[indices[indptr[i]:indptr[i+1]]]*x[indices[indptr[i]:indptr[i+1]]])
213        tmp = 0.0
214        for j in range(indptr[i], indptr[i+1]):
215            ind = indices[j]
216            tmp += data[j]*x[ind]
217        b[i] = tmp
218       
219    return b
220"""
221   
222class MasterAlgOptions(OptionBase):
223    """
224    Options for solving coupled FMI 2 CS FMUs.
225
226    Options::
227
228        step_size --
229            Specfies the global step-size to be used for simulating
230            the coupled system.
231            Default: 0.01
232           
233        initialize --
234            If set to True, the initializing algorithm defined in the FMU models
235            are invoked, otherwise it is assumed the user have manually initialized
236            all models.
237            Default is True.
238           
239        block_initialization --
240            If set to True, the initialization algorithm computes the evaluation
241            order of the FMUs and tries to resolve algebraic loops by this
242            evaluation order.
243            Default is False.
244       
245        extrapolation_order --
246            Defines the extrapolation used in the simulation.
247            Default is 0 (constant extrapolation).
248       
249        smooth_coupling --
250            Defines if the extrapolation should be smoothen, i.e. the input
251            values are adapted so that they are C^0 instead of C^(-1) in case
252            extrapolation_order is > 0.
253            Default is True
254       
255        linear_correction --
256            Defines if linear correction should be used during the simulation.
257            Note that this increases the simulation robustness in case of
258            algebraic loops.
259            Default is True
260
261        execution --
262            Defines if the models are to be evaluated in parallel (note that it
263            is not an algorithm change, just an evaluation execution within
264            the same algorithm). Note that it requires that PyFMI has been
265            installed with OpenMP.
266            Default is serial
267           
268        num_threads --
269            Defines the number of threads used when the execution is set
270            to parallel.
271            Default: Number of cores / OpenMP environment variable
272           
273        error_controlled --
274            Defines if the algorithm should adapt the step-size during
275            the simulation. Note requires that all FMUs support save/get
276            state.
277            Default: False
278           
279        atol --
280            Defines the absolute tolerance used when an error controlled
281            simulation is performed.
282            Default: 1e-4
283
284        rtol --
285            Defines the relative tolerance used when an error controlled
286            simulation is performed.
287            Default: 1e-4
288           
289        maxh --
290            Defines the maximum step-size allowed to be used together
291            with an error controlled simulation.
292            Default: 0.0 (i.e. inactive)
293
294        result_file_name --
295            Specifies the name of the file where the simulation result is
296            written. Setting this option to an empty string results in a default
297            file name that is based on the name of the model class.
298            Default: Empty string
299
300        result_handling --
301            Specifies how the result should be handled. Either stored to
302            file or stored in memory. One can also use a custom handler.
303            Available options: "file", "memory", "csv", "custom"
304            Default: "file"
305
306        result_handler --
307            The handler for the result. Depending on the option in
308            result_handling this either defaults to ResultHandlerFile
309            or ResultHandlerMemory. If result_handling custom is choosen
310            This MUST be provided.
311            Default: None
312
313        filter --
314            A filter for choosing which variables to actually store
315            result for. The syntax can be found in
316            http://en.wikipedia.org/wiki/Glob_%28programming%29 . An
317            example is filter = "*der" , stor all variables ending with
318            'der'. Can also be a list. Note that there should be one
319            filter for each model.
320            Default: None
321
322
323    """
324    def __init__(self, master, *args, **kw):
325        _defaults= {
326        "initialize" : True,
327        "local_rtol" : 1e-6,
328        "rtol"       : 1e-4,
329        "atol"       : 1e-4,
330        "step_size"  : 0.01,
331        "maxh"       : 0.0,
332        "filter"     : dict((model,None) for model in master.models),
333        "result_handling"     : "file",
334        "result_handler"      : None,
335        "inputs"              : None,
336        "linear_correction"   : True,
337        "error_controlled"    : False if master.support_storing_fmu_states else False,
338        "logging"             : False,
339        "extrapolation_order" : 0, #Constant
340        "store_step_before_update" : False,
341        "smooth_coupling"             : True,
342        "execution"                : "serial",
343        "block_initialization"     : False,
344        "block_initialization_type" : "greedy",
345        "block_initialization_order" : None,
346        "experimental_output_derivative": False,
347        "experimental_finite_difference_D": False,
348        "experimental_output_solve":False,
349        "force_finite_difference_outputs": False,
350        "num_threads":None}
351        super(MasterAlgOptions,self).__init__(_defaults)
352        self._update_keep_dict_defaults(*args, **kw)
353
354cdef class Master:
355    cdef public list connections, models
356    cdef public dict statistics, models_id_mapping
357    cdef public object opts
358    cdef public object models_dict, L
359    cdef public object I
360    cdef public object y_prev, yd_prev, input_traj
361    cdef public object DL_prev
362    cdef public int algebraic_loops, storing_fmu_state
363    cdef public int error_controlled, linear_correction
364    cdef public int _support_directional_derivatives, _support_storing_fmu_states, _support_interpolate_inputs, _max_output_derivative_order
365    cdef double rtol, atol, current_step_size
366    cdef public object y_m1, yd_m1, u_m1, ud_m1, udd_m1
367    cdef FMIL.fmi2_import_t** fmu_adresses
368    cdef public int _len_inputs, _len_outputs
369    cdef public int _len_derivatives
370    cdef public list _storedDrow, _storedDcol
371    cdef public np.ndarray _array_one
372    cdef public object _D
373    cdef public dict elapsed_time
374    cdef public dict elapsed_time_init
375    cdef public dict _error_data
376    cdef public int _display_counter
377    cdef public object _display_progress
378    cdef public double _time_integration_start
379   
380    def __init__(self, models, connections):
381        """
382        Initializes the master algorithm.
383       
384        Parameters::
385       
386            models 
387                    - A list of models that are to be simulated.
388                      Needs to be a subclass of FMUModelCS.
389                     
390            connection 
391                    - Specifices the connection between the models.
392                       
393                    - model_begin.variable -> model_accept.variable
394                      [(model_source,"beta",model_destination,"y"),(...)]
395       
396        """
397        if not isinstance(models, list):
398            raise fmi.FMUException("The models should be provided as a list.")
399        for model in models:
400            if not isinstance(model, fmi.FMUModelCS2):
401                raise fmi.FMUException("The Master algorithm currently only supports CS 2.0 FMUs.")
402        self.fmu_adresses = <FMIL.fmi2_import_t**>FMIL.malloc(len(models)*sizeof(FMIL.fmi2_import_t*))
403       
404        self.connections = connections
405        self.models = models
406        self.models_dict = OrderedDict((model,{"model": model, "result": None, "external_input": None, 
407                                               "local_input": [], "local_input_vref": [], "local_input_len": 0,
408                                               "local_state": [], "local_state_vref": [],
409                                               "local_derivative": [], "local_derivative_vref": [],
410                                               "local_output": [], "local_output_vref": [], "local_output_len": 0,
411                                               "local_output_range_array": None,
412                                               "direct_dependence": []}) for model in models)
413        self.models_id_mapping = {str(id(model)): model for model in models}
414        self.elapsed_time = {model: 0.0 for model in models}
415        self.elapsed_time_init = {model: 0.0 for model in models}
416        self.elapsed_time["result_handling"] = 0.0
417        self._display_counter = 1
418        self._display_progress = True
419       
420        self.statistics = {}
421        self.statistics["nsteps"] = 0
422        self.statistics["nreject"] = 0
423       
424        #Initialize internal variables
425        self._support_directional_derivatives = -1
426        self._support_storing_fmu_states = -1
427        self._support_interpolate_inputs = -1
428        self._max_output_derivative_order = -1
429        self._len_inputs = 0
430        self._len_outputs = 0
431        self._len_derivatives = 0
432        self._array_one = np.array([1.0])
433        self._D = None
434       
435        self.error_controlled = 0
436        self.linear_correction = 1
437       
438        self.check_support_storing_fmu_state()
439       
440        self.connection_setup(connections)
441        self.verify_connection_variables()
442        self.check_algebraic_loops()
443        self.set_model_order()
444        self.define_connection_matrix()
445       
446        self.y_prev = None
447        self.input_traj = None
448        self.I = sp.eye(self._len_inputs, self._len_outputs, format="csr") #y = Cx + Du , u = Ly -> DLy   DL[inputsXoutputs]
449       
450        self._error_data = {"time":[], "error":[], "step-size":[], "rejected":[]}
451   
452    def __del__(self):
453        FMIL.free(self.fmu_adresses)
454   
455    cdef set_last_y(self, np.ndarray y):
456        self.y_m1 = y.copy()
457    cdef get_last_y(self):
458        return self.y_m1
459    cdef set_last_yd(self, np.ndarray yd):
460        self.yd_m1 = yd.copy() if yd is not None else None
461    cdef get_last_yd(self):
462        return self.yd_m1
463    cdef set_last_us(self, np.ndarray u, np.ndarray ud=None, np.ndarray udd=None):
464        self.u_m1 = u.copy()
465        self.ud_m1 = ud.copy() if ud is not None else None
466        self.udd_m1 = udd.copy() if udd is not None else None
467    cdef get_last_us(self):
468        return self.u_m1, self.ud_m1, self.udd_m1
469    cdef set_current_step_size(self, double step_size):
470        self.current_step_size = step_size
471    cdef double get_current_step_size(self):
472        return self.current_step_size
473       
474    def report_solution(self, double cur_time):
475       
476        store_communication_point(self.models_dict)
477       
478        if self._display_progress:
479            if ( timer() - self._time_integration_start) > self._display_counter*10:
480                self._display_counter += 1
481               
482                sys.stdout.write(" Simulation time: %e" % cur_time)
483                sys.stdout.write('\r')
484                sys.stdout.flush()
485       
486    def set_model_order(self):
487        i = 0
488        for model in self.models_dict.keys():
489            self.models_dict[model]["order"] = i
490            i = i+1
491        for model in self.models_dict.keys():
492            self.models[self.models_dict[model]["order"]] = model
493           
494    def copy_fmu_addresses(self):
495        for model in self.models_dict.keys():
496            self.fmu_adresses[self.models_dict[model]["order"]] = (<FMUModelCS2>model)._fmu
497           
498    def define_connection_matrix(self):
499        cdef list data = []
500        cdef list row = []
501        cdef list col = []
502
503        start_index_inputs      = 0
504        start_index_outputs     = 0     
505        start_index_states      = 0
506        start_index_derivatives = 0
507        for model in self.models_dict.keys():
508            self.models_dict[model]["global_index_inputs"]      = start_index_inputs
509            self.models_dict[model]["global_index_outputs"]     = start_index_outputs
510            self.models_dict[model]["global_index_states"]      = start_index_states
511            self.models_dict[model]["global_index_derivatives"] = start_index_derivatives
512           
513            start_index_inputs      += len(self.models_dict[model]["local_input"])
514            start_index_outputs     += len(self.models_dict[model]["local_output"])
515            start_index_states      += len(self.models_dict[model]["local_state"])
516            start_index_derivatives += len(self.models_dict[model]["local_derivative"])
517       
518        for connection in self.connections:
519            src = connection[0]; src_var = connection[1]
520            dst = connection[2]; dst_var = connection[3]
521            data.append(1)
522            row.append(self.models_dict[dst]["global_index_inputs"]+self.models_dict[dst]["local_input"].index(dst_var))
523            col.append(self.models_dict[src]["global_index_outputs"]+self.models_dict[src]["local_output"].index(src_var))
524           
525        self.L = sp.csr_matrix((data, (row, col)), (len(self.connections),len(self.connections)), dtype=np.float64)
526       
527    cpdef compute_global_D(self):
528        cdef list data = []
529        cdef list row = []
530        cdef list col = []
531        cdef int i, nlocal, status
532
533        for model in self.models_dict.keys():
534            nlocal = self.models_dict[model]["local_input_len"]
535            #v = [0.0]*nlocal
536            for i in range(nlocal):
537                #local_D = model.get_directional_derivative([self.models_dict[model]["local_input_vref"][i]],self.models_dict[model]["local_output_vref"], [1.0])
538                local_D = np.empty(self.models_dict[model]["local_output_len"])
539                #status = (<FMUModelCS2>model)._get_directional_derivative(np.array([self.models_dict[model]["local_input_vref"][i]]),self.models_dict[model]["local_output_vref_array"], self._array_one, local_D)
540                if self.opts["experimental_finite_difference_D"]:
541                    up = (<FMUModelCS2>model).get_real(self.models_dict[model]["local_input_vref_array"][i:i+1])
542                    eps = max(abs(up), 1.0)
543                    yp = (<FMUModelCS2>model).get_real(self.models_dict[model]["local_output_vref_array"])
544                    (<FMUModelCS2>model).set_real(self.models_dict[model]["local_input_vref_array"][i:i+1], up+eps)
545                    local_D = ((<FMUModelCS2>model).get_real(self.models_dict[model]["local_output_vref_array"]) - yp)/eps
546                    (<FMUModelCS2>model).set_real(self.models_dict[model]["local_input_vref_array"][i:i+1], up)
547                else:
548                    status = (<FMUModelCS2>model)._get_directional_derivative(self.models_dict[model]["local_input_vref_array"][i:i+1],self.models_dict[model]["local_output_vref_array"], self._array_one, local_D)
549           
550                    if status != 0: raise fmi.FMUException("Failed to get the directional derivatives while computing the global D matrix.")
551                data.extend(local_D)
552               
553                if self._storedDrow is None and self._storedDcol is None:
554                    col.extend([self.models_dict[model]["global_index_inputs"]+i]*len(local_D))
555                    #row.extend(np.array([self.models_dict[model]["global_index_outputs"]]*self.models_dict[model]["local_output_len"])+np.array(range(self.models_dict[model]["local_output_len"])))
556                    row.extend(np.array([self.models_dict[model]["global_index_outputs"]]*self.models_dict[model]["local_output_len"])+self.models_dict[model]["local_output_range_array"])
557       
558        if self._storedDrow is None and self._storedDcol is None:
559            self._storedDrow = row
560            self._storedDcol = col
561        else:
562            row = self._storedDrow
563            col = self._storedDcol
564       
565        if self._D is None:
566            self._D = sp.csr_matrix((data, (row, col)))#, (len(col),len(row)))
567        else:
568            self._D.data = np.array(data, dtype=np.float64)
569           
570        return self._D
571           
572   
573    def compute_global_C(self):
574        cdef list data = []
575        cdef list row = []
576        cdef list col = []
577
578        for model in self.models_dict.keys():
579            if model.get_generation_tool() != "JModelica.org":
580                return None
581            v = [0.0]*len(self.models_dict[model]["local_state_vref"])
582            for i in range(len(v)):
583                local_C = model.get_directional_derivative([self.models_dict[model]["local_state_vref"][i]],self.models_dict[model]["local_output_vref"], [1.0])
584                data.extend(local_C)
585                col.extend([self.models_dict[model]["global_index_states"]+i]*len(local_C))
586                row.extend(np.array([self.models_dict[model]["global_index_outputs"]]*len(self.models_dict[model]["local_output_vref"]))+np.array(range(len(self.models_dict[model]["local_output_vref"]))))
587        return sp.csr_matrix((data, (row, col)))
588       
589    def compute_global_A(self):
590        cdef list data = []
591        cdef list row = []
592        cdef list col = []
593
594        for model in self.models_dict.keys():
595            if model.get_generation_tool() != "JModelica.org":
596                return None
597            v = [0.0]*len(self.models_dict[model]["local_state_vref"])
598            for i in range(len(v)):
599                local_A = model.get_directional_derivative([self.models_dict[model]["local_state_vref"][i]],self.models_dict[model]["local_derivative_vref"], [1.0])
600                data.extend(local_A)
601                col.extend([self.models_dict[model]["global_index_states"]+i]*len(local_A))
602                row.extend(np.array([self.models_dict[model]["global_index_derivatives"]]*len(self.models_dict[model]["local_derivative_vref"]))+np.array(range(len(self.models_dict[model]["local_derivative_vref"]))))
603        return sp.csr_matrix((data, (row, col)))
604       
605    def compute_global_B(self):
606        cdef list data = []
607        cdef list row = []
608        cdef list col = []
609
610        for model in self.models_dict.keys():
611            if model.get_generation_tool() != "JModelica.org":
612                return None
613            v = [0.0]*len(self.models_dict[model]["local_input_vref"])
614            for i in range(len(v)):
615                local_B = model.get_directional_derivative([self.models_dict[model]["local_input_vref"][i]],self.models_dict[model]["local_derivative_vref"], [1.0])
616                data.extend(local_B)
617                col.extend([self.models_dict[model]["global_index_inputs"]+i]*len(local_B))
618                row.extend(np.array([self.models_dict[model]["global_index_derivatives"]]*len(self.models_dict[model]["local_derivative_vref"]))+np.array(range(len(self.models_dict[model]["local_derivative_vref"]))))
619        return sp.csr_matrix((data, (row, col)))
620       
621    def connection_setup(self, connections):
622        for connection in connections:
623            self.models_dict[connection[0]]["local_output"].append(connection[1])
624            self.models_dict[connection[0]]["local_output_vref"].append(connection[0].get_variable_valueref(connection[1]))
625            self.models_dict[connection[2]]["local_input"].append(connection[3])
626            self.models_dict[connection[2]]["local_input_vref"].append(connection[2].get_variable_valueref(connection[3]))
627           
628        for model in self.models_dict.keys():
629            self.models_dict[model]["local_input_len"] = len(self.models_dict[model]["local_input"])
630            self.models_dict[model]["local_output_len"] = len(self.models_dict[model]["local_output"])
631            self.models_dict[model]["local_output_range_array"] = np.array(range(self.models_dict[model]["local_output_len"]))
632            self.models_dict[model]["local_output_vref_array"] = np.array(self.models_dict[model]["local_output_vref"], dtype=np.uint32)
633            self.models_dict[model]["local_input_vref_array"] = np.array(self.models_dict[model]["local_input_vref"], dtype=np.uint32)
634            self.models_dict[model]["local_input_vref_ones"] = np.ones(self.models_dict[model]["local_input_len"], dtype=np.int32)
635            self.models_dict[model]["local_input_vref_twos"] = 2*np.ones(self.models_dict[model]["local_input_len"], dtype=np.int32)
636            self.models_dict[model]["local_output_vref_ones"] = np.ones(self.models_dict[model]["local_output_len"], dtype=np.int32)
637            self._len_inputs  += self.models_dict[model]["local_input_len"]
638            self._len_outputs += self.models_dict[model]["local_output_len"]
639           
640            if model.get_generation_tool() == "JModelica.org":
641                self.models_dict[model]["local_state"]           = model.get_states_list().keys()
642                self.models_dict[model]["local_state_vref"]      = [var.value_reference for var in model.get_states_list().values()]
643                self.models_dict[model]["local_derivative"]      = model.get_derivatives_list().keys()
644                self.models_dict[model]["local_derivative_vref"] = [var.value_reference for var in model.get_derivatives_list().values()]
645                self.models_dict[model]["local_derivative_vref_array"] = np.array(self.models_dict[model]["local_derivative_vref"], dtype=np.uint32)
646                self.models_dict[model]["local_derivative_len"] = len(self.models_dict[model]["local_derivative"])
647                self._len_derivatives += self.models_dict[model]["local_derivative_len"]
648               
649    def verify_connection_variables(self):
650        for model in self.models_dict.keys():
651            for output in self.models_dict[model]["local_output"]:
652                if model.get_variable_causality(output) != fmi.FMI2_OUTPUT:
653                    raise fmi.FMUException("The connection variable " + output + " in model " + model.get_name() + " is not an output. ")
654            for input in self.models_dict[model]["local_input"]:
655                if model.get_variable_causality(input) != fmi.FMI2_INPUT:
656                    raise fmi.FMUException("The connection variable " + input + " in model " + model.get_name() + " is not an input. ")
657                   
658    def check_algebraic_loops(self):
659        """
660        Simplified check for algebraic loops in simulation mode due to
661        the limited capacity of solving the loops
662        """
663        self.algebraic_loops = 0
664       
665        for model in self.models_dict.keys():
666            output_state_dep, output_input_dep = model.get_output_dependencies()
667            for local_output in self.models_dict[model]["local_output"]:
668                for local_input in self.models_dict[model]["local_input"]:
669                    if local_input in output_input_dep[local_output]:
670                        self.models_dict[model]["direct_dependence"].append((local_input, local_output))
671                        self.algebraic_loops = 1
672                        #break
673                if self.algebraic_loops:
674                    pass
675                    #break
676            if self.algebraic_loops:
677                pass
678                #break
679                   
680        if self.algebraic_loops:
681            for model in self.models_dict.keys():
682                if model.get_capability_flags()["providesDirectionalDerivatives"] is False:
683                    warnings.warn("The model, " + model.get_name() + ", does not support " 
684                                "directional derivatives which is necessary in-case of an algebraic loop. The simulation might become unstable...")
685       
686        return self.algebraic_loops
687       
688    def check_support_storing_fmu_state(self):
689        self.storing_fmu_state = 1
690        for model in self.models_dict.keys():
691            if model.get_capability_flags()["canGetAndSetFMUstate"] is False:
692                self.storing_fmu_state= 0
693                break
694        return self.storing_fmu_state
695       
696    cpdef np.ndarray get_connection_outputs(self):
697        cdef int i = 0, inext = 0
698        cdef np.ndarray y = np.empty((self._len_outputs))
699        #for model in self.models_dict.keys():
700        for model in self.models:
701            #y.extend(model.get(self.models_dict[model]["local_output"]))
702            #y.extend(model.get_real(self.models_dict[model]["local_output_vref"]))
703            i = self.models_dict[model]["global_index_outputs"]
704            inext = i + self.models_dict[model]["local_output_len"]
705            y[i:inext] = (<FMUModelCS2>model).get_real(self.models_dict[model]["local_output_vref_array"])
706            i = inext
707           
708        #return np.array(y)
709        return y.reshape(-1,1)
710       
711    cpdef np.ndarray _get_derivatives(self):
712        cdef int i = 0, inext = 0
713        cdef np.ndarray xd = np.empty((self._len_derivatives))
714       
715        for model in self.models_dict.keys():
716            if model.get_generation_tool() != "JModelica.org":
717                return None
718
719        for model in self.models:
720            i = self.models_dict[model]["global_index_derivatives"]
721            inext = i + self.models_dict[model]["local_derivative_len"]
722            xd[i:inext] = (<FMUModelCS2>model).get_real(self.models_dict[model]["local_derivative_vref_array"])
723            i = inext
724
725        return xd.reshape(-1,1)
726       
727    cpdef np.ndarray get_specific_connection_outputs(self, model, np.ndarray mask, np.ndarray yout):
728        cdef int j = 0
729        cdef np.ndarray ytmp = (<FMUModelCS2>model).get_real(self.models_dict[model]["local_output_vref_array"][mask])
730        for i, flag in enumerate(mask):
731            if flag:
732                yout[i+self.models_dict[model]["global_index_outputs"]] = ytmp[j]
733                j = j + 1
734       
735    cpdef get_connection_derivatives(self, np.ndarray y_cur):
736        #cdef list yd = []
737        cdef int i = 0, inext = 0, status = 0
738        cdef np.ndarray[FMIL.fmi2_real_t, ndim=1, mode='c']  yd    = np.empty((self._len_outputs))
739        cdef np.ndarray[FMIL.fmi2_real_t, ndim=1, mode='c']  ydtmp = np.empty((self._len_outputs))
740        cdef np.ndarray y_last = None
741       
742        if self.opts["extrapolation_order"] > 0:
743            if self.max_output_derivative_order > 0 and not self.opts["force_finite_difference_outputs"]:
744                for model in self.models_dict.keys():
745                    #yd.extend(model.get_output_derivatives(self.models_dict[model]["local_output"], 1))
746                    inext = i + self.models_dict[model]["local_output_len"]
747                    status = (<FMUModelCS2>model)._get_output_derivatives(self.models_dict[model]["local_output_vref_array"], ydtmp, self.models_dict[model]["local_output_vref_ones"])
748                    if status != 0: raise fmi.FMUException("Failed to get the output derivatives.")
749                    yd[i:inext] = ydtmp[:inext-i]
750                    i = inext
751               
752                return self.correct_output_derivative(yd.reshape(-1,1))
753                #return yd.reshape(-1,1)
754               
755            else:
756                       
757                if self.opts["experimental_output_derivative"]:
758                    JM_FMUS = True
759                    for model in self.models_dict.keys():
760                        if model.get_generation_tool() != "JModelica.org":
761                            JM_FMUS = False
762                            break
763                    if JM_FMUS:
764                        C = self.compute_global_C()
765                        D = self.compute_global_D()
766                        u,ud,udd = self.get_last_us()
767                        xd = self._get_derivatives()
768                        if ud is not None:
769                            if udd is not None:
770                                return C.dot(xd)+D.dot(ud+self.get_current_step_size()*udd)
771                            else:
772                                return C.dot(xd)+D.dot(ud)
773                        else: #First step
774                            return splin.spsolve((self.I-D.dot(self.L)),C.dot(xd)).reshape((-1,1))
775
776                y_last = self.get_last_y()
777                if y_last is not None:
778                    return (y_cur - y_last)/self.get_current_step_size()
779                else:
780                    return None
781        else:
782            return None
783           
784    cpdef get_connection_second_derivatives(self, np.ndarray yd_cur):
785
786        cdef int i = 0, inext = 0, status = 0
787        cdef np.ndarray[FMIL.fmi2_real_t, ndim=1, mode='c']  ydd    = np.empty((self._len_outputs))
788        cdef np.ndarray[FMIL.fmi2_real_t, ndim=1, mode='c']  yddtmp = np.empty((self._len_outputs))
789        cdef np.ndarray yd_last = None
790       
791        if self.opts["extrapolation_order"] > 1:
792            if self.max_output_derivative_order > 1 and not self.opts["force_finite_difference_outputs"]:
793                for model in self.models_dict.keys():
794                    inext = i + self.models_dict[model]["local_output_len"]
795                    status = (<FMUModelCS2>model)._get_output_derivatives(self.models_dict[model]["local_output_vref_array"], yddtmp, self.models_dict[model]["local_output_vref_twos"])
796                    if status != 0: raise fmi.FMUException("Failed to get the output derivatives of second order.")
797                    ydd[i:inext] = yddtmp[:inext-i]
798                    i = inext
799               
800                return self.correct_output_second_derivative(ydd.reshape(-1,1))
801               
802            else:
803                if self.opts["experimental_output_derivative"]:
804                    JM_FMUS = True
805                    for model in self.models_dict.keys():
806                        if model.get_generation_tool() != "JModelica.org":
807                            JM_FMUS = False
808                            break
809                    if JM_FMUS:
810                        A = self.compute_global_A()
811                        B = self.compute_global_B()
812                        C = self.compute_global_C()
813                        D = self.compute_global_D()
814                        u,ud,udd = self.get_last_us()
815                        xd = self._get_derivatives()
816                        if ud is not None and udd is not None:
817                            return C.dot(A.dot(xd))+C.dot(B.dot(ud+self.get_current_step_size()*udd))+D.dot(udd)
818                        else: #First step
819                            return splin.spsolve((self.I-D.dot(self.L)),C.dot(A.dot(xd)+B.dot(self.L.dot(yd_cur)))).reshape((-1,1))
820               
821                yd_last = self.get_last_yd()
822                if yd_last is not None:
823                    return (yd_cur - yd_last)/self.get_current_step_size()
824                else:
825                    return None
826        else:
827            return None
828           
829    cpdef set_connection_inputs(self, np.ndarray u, np.ndarray ud=None, np.ndarray udd=None):
830        cdef int i = 0, inext, status
831        #for model in self.models_dict.keys():
832        u = u.ravel()
833        for model in self.models:
834            i = self.models_dict[model]["global_index_inputs"] #MIGHT BE WRONG
835            inext = i + self.models_dict[model]["local_input_len"]
836            #model.set(self.models_dict[model]["local_input"], u[i:inext])
837            (<FMUModelCS2>model).set_real(self.models_dict[model]["local_input_vref_array"], u[i:inext])
838           
839            if ud is not None: #Set the input derivatives
840                ud = ud.ravel()
841                #model.set_input_derivatives(self.models_dict[model]["local_input"], ud[i:inext], 1)
842                status = (<FMUModelCS2>model)._set_input_derivatives(self.models_dict[model]["local_input_vref_array"], ud[i:inext], self.models_dict[model]["local_input_vref_ones"])
843                if status != 0: raise fmi.FMUException("Failed to set the first order input derivatives.")
844               
845            if udd is not None: #Set the input derivatives
846                udd = udd.ravel()
847                status = (<FMUModelCS2>model)._set_input_derivatives(self.models_dict[model]["local_input_vref_array"], udd[i:inext], self.models_dict[model]["local_input_vref_twos"])
848                if status != 0: raise fmi.FMUException("Failed to set the second order input derivatives.")
849           
850            i = inext
851           
852    cpdef set_specific_connection_inputs(self, model, np.ndarray mask, np.ndarray u):
853        cdef int i = self.models_dict[model]["global_index_inputs"]
854        cdef int inext = i + self.models_dict[model]["local_input_len"]
855        cdef np.ndarray usliced = u.ravel()[i:inext]
856        (<FMUModelCS2>model).set_real(self.models_dict[model]["local_input_vref_array"][mask], usliced[mask])
857
858    cpdef correct_output_second_derivative(self, np.ndarray ydd):
859        if self.linear_correction and self.algebraic_loops and self.support_directional_derivatives:
860            raise NotImplementedError
861            """
862            D = self.compute_global_D()
863            DL = D.dot(self.L)
864           
865            if self.opts["extrapolation_order"] > 0:
866                uold, udold = self.get_last_us()
867                uhat = udold if udold is not None else np.zeros(np.array(yd).shape)
868               
869                z = yd - D.dot(uhat)
870           
871            yd = splin.spsolve((self.I-DL),z).reshape((-1,1))
872            """
873        return ydd
874
875    cpdef correct_output_derivative(self, np.ndarray yd):
876        if self.linear_correction and self.algebraic_loops and self.support_directional_derivatives:
877            D = self.compute_global_D()
878            DL = D.dot(self.L)
879           
880            if self.opts["extrapolation_order"] > 0:
881                uold, udold, uddold = self.get_last_us()
882                uhat = udold if udold is not None else np.zeros(np.array(yd).shape)
883               
884                z = yd - D.dot(uhat)
885           
886            yd = splin.spsolve((self.I-DL),z).reshape((-1,1))
887
888        return yd
889   
890    cpdef correct_output(self, np.ndarray y, y_prev=None):
891        if self.algebraic_loops and self.opts["experimental_output_solve"]:
892            JM_FMUS = True
893            for model in self.models_dict.keys():
894                if model.get_generation_tool() != "JModelica.org":
895                    JM_FMUS = False
896                    break
897            if JM_FMUS:
898                if USE_ROOT:
899                    res = sopt.root(init_f, y, args=(self))
900                else:
901                    res = sopt.fsolve(init_f, y, args=(self))
902                if not res["success"]:
903                    print res
904                    raise Exception("Failed to converge the output system.")
905                return res["x"].reshape(-1,1)
906               
907        if self.linear_correction and self.algebraic_loops and y_prev is not None:# and self.support_directional_derivatives:
908            D = self.compute_global_D()
909            DL = D.dot(self.L)
910           
911            if self.opts["extrapolation_order"] > 0:
912                uold, udold, uddold = self.get_last_us()
913                uhat = uold + (self.get_current_step_size()*udold if udold is not None else 0.0)
914               
915                z = y - D.dot(uhat)
916                #z = y - matvec(D,uhat.ravel())
917            else:
918               
919                z = y - DL.dot(y_prev)
920                #z = y - matvec(DL, y_prev.ravel())
921           
922            y = splin.spsolve((self.I-DL),z).reshape((-1,1))
923            #y = splin.lsqr((sp.eye(*DL.shape)-DL),z)[0].reshape((-1,1))
924
925        elif self.algebraic_loops and self.support_directional_derivatives:
926            pass
927           
928        return y
929   
930    cpdef modify_input(self, np.ndarray y, np.ndarray yd, np.ndarray ydd):
931        cdef double h = self.get_current_step_size()
932       
933        u    = self.L.dot(y)
934        #u    = matvec(self.L,y.ravel())
935        ud   = self.L.dot(yd) if yd is not None else None
936        #ud   = matvec(self.L,yd.ravel()) if yd is not None else None
937        udd  = self.L.dot(ydd) if ydd is not None else None
938       
939        if self.opts["extrapolation_order"] > 0 and self.opts["smooth_coupling"]:
940            uold, udold, uddold = self.get_last_us()
941            uhat = uold + (h*udold if udold is not None else 0.0)
942       
943            udhat = (u-uhat)/h+ud if ud is not None else ud
944           
945            u = uhat
946            ud = udhat
947
948        return u, ud, udd
949       
950   
951    cpdef exchange_connection_data(self):
952        #u = Ly
953        cdef np.ndarray y = self.get_connection_outputs()
954        cdef np.ndarray u
955       
956        y   = self.correct_output(y, self.y_prev)
957        yd  = self.get_connection_derivatives(y)
958        ydd = self.get_connection_second_derivatives(yd)
959       
960        self.y_prev = y.copy()
961        self.yd_prev = yd.copy() if yd is not None else None
962       
963        u, ud, udd = self.modify_input(y, yd, ydd)
964       
965        self.set_connection_inputs(u, ud=ud, udd=udd)
966       
967        self.set_last_us(u, ud, udd)
968       
969        return y, yd, u
970   
971    def initialize(self, double start_time, double final_time, object opts):
972        self.set_input(start_time)
973       
974        if opts["block_initialization"]:
975            order, blocks, compressed = self.compute_evaluation_order(opts["block_initialization_type"], order=opts["block_initialization_order"])
976            model_in_init_mode = {model:False for model in self.models}
977           
978            #Global outputs vector
979            y = np.zeros((self._len_outputs))
980           
981            for block in blocks:
982                if len(block["inputs"]) == 0: #No inputs in this block, only outputs
983                    for model in block["outputs"].keys():
984                        #If the model has not entered initialization mode, enter
985                        if model_in_init_mode[model] is False:
986                            enter_initialization_mode([model], start_time, final_time, opts, self.elapsed_time_init)
987                            model_in_init_mode[model] = True
988                       
989                        #Get the outputs
990                        time_start = timer()
991                        self.get_specific_connection_outputs(model, block["outputs_mask"][model], y)
992                        self.elapsed_time_init[model] += timer() - time_start
993                       
994                elif len(block["outputs"]) == 0: #No outputs in this block
995                   
996                    #Compute current global input vector
997                    u = self.L.dot(y.reshape(-1,1))
998                    for model in block["inputs"].keys(): #Set the inputs
999                        #print "Model: ", model
1000                        #print "Inputs: ", block["inputs"]
1001                        #print "Mask: ", block["inputs_mask"]
1002                        self.set_specific_connection_inputs(model, block["inputs_mask"][model], u)
1003                       
1004                else: #Both (algebraic loop)
1005                   
1006                    #Assert models has entered initialization mode
1007                    for model in block["inputs"].keys()+block["outputs"].keys(): #Possible only need outputs?
1008                        if model_in_init_mode[model] is False:
1009                            enter_initialization_mode([model], start_time, final_time, opts, self.elapsed_time_init)
1010                            model_in_init_mode[model] = True
1011                           
1012                    #Get start values
1013                    for model in block["outputs"].keys():
1014                        self.get_specific_connection_outputs(model, block["outputs_mask"][model], y)
1015                   
1016                    if USE_ROOT:
1017                        res = sopt.root(init_f_block, y[block["global_outputs_mask"]], args=(self,block))
1018                    else:
1019                        res = sopt.fsolve(init_f_block, y[block["global_outputs_mask"]], args=(self,block))
1020                    if not res["success"]:
1021                        print res
1022                        raise Exception("Failed to converge the initialization system.")
1023                   
1024                    y[block["global_outputs_mask"]] = res["x"]
1025                    u = self.L.dot(y.reshape(-1,1))
1026                   
1027                    for model in block["inputs"].keys():
1028                        #Set the inputs
1029                        self.set_specific_connection_inputs(model, block["inputs_mask"][model], u)
1030                   
1031            #Assert that all models has entered initialization mode       
1032            for model in self.models:
1033                if model_in_init_mode[model] is False:
1034                    enter_initialization_mode([model], start_time, final_time, opts, self.elapsed_time_init)
1035                    model_in_init_mode[model] = True
1036        else:
1037            enter_initialization_mode(self.models, start_time, final_time, opts, self.elapsed_time_init)
1038           
1039            if self.algebraic_loops: #If there is an algebraic loop, solve the resulting system
1040                if self.support_directional_derivatives:
1041                    if USE_ROOT:
1042                        res = sopt.root(init_f, self.get_connection_outputs(), args=(self), jac=init_jac)
1043                    else:
1044                        res = sopt.fsolve(init_f, self.get_connection_outputs(), args=(self), jac=init_jac)
1045                else:
1046                    if USE_ROOT:
1047                        res = sopt.root(init_f, self.get_connection_outputs(), args=(self))
1048                    else:
1049                        res = sopt.fsolve(init_f, self.get_connection_outputs(), args=(self))
1050                if not res["success"]:
1051                    print res
1052                    raise Exception("Failed to converge the initialization system.")
1053                u = self.L.dot(res["x"].reshape(-1,1))
1054                self.set_connection_inputs(u)
1055            else:
1056                y = self.get_connection_outputs()
1057                u = self.L.dot(y)
1058                self.set_connection_inputs(u)
1059       
1060        exit_initialization_mode(self.models, self.elapsed_time_init)
1061       
1062        #Store the outputs
1063        self.y_prev = self.get_connection_outputs().copy()
1064        self.set_last_y(self.y_prev)
1065       
1066    def initialize_result_objects(self, opts):
1067        i = 0
1068        for model in self.models_dict.keys():
1069            if opts["result_handling"] == "file":
1070                result_object = ResultHandlerFile(model)
1071            elif opts["result_handling"] == "none":
1072                result_object = ResultHandlerDummy(model)
1073            else:
1074                raise fmi.FMUException("Currently only writing result to file and none is supported.")
1075            local_opts = FMICSAlgOptions()
1076            local_opts["result_file_name"] = model.get_identifier()+'_'+str(i)+'_result.txt'
1077            local_opts["filter"] = opts["filter"][model]
1078           
1079            result_object.set_options(local_opts)
1080            result_object.simulation_start()
1081            result_object.initialize_complete()
1082           
1083            i = i + 1
1084           
1085            self.models_dict[model]["result"] = result_object
1086           
1087    def jacobi_algorithm(self, double start_time, double final_time, object opts):
1088        cdef double step_size = opts["step_size"]
1089        cdef int calling_setting = SERIAL if opts["execution"] != "parallel" else PARALLEL
1090        cdef double tcur, step_size_old
1091        cdef double error
1092        cdef dict states = None
1093        cdef np.ndarray ycur, ucur
1094       
1095        if self.error_controlled:
1096            tcur = start_time
1097           
1098            y_old = self.get_connection_outputs()
1099           
1100            while tcur < final_time:
1101                #Store FMU states
1102                states = get_fmu_states(self.models, states)
1103               
1104                #Set input
1105                self.set_input(tcur)
1106                #Take a full step
1107                perform_do_step(self.models, self.elapsed_time, self.fmu_adresses, tcur, step_size, False, calling_setting)
1108               
1109                y_full = self.correct_output(self.get_connection_outputs(), y_old)
1110               
1111                #Restore FMU states
1112                set_fmu_states(states)
1113               
1114                #Set input
1115                self.set_input(tcur)
1116                #Take a half step
1117                perform_do_step(self.models, self.elapsed_time, self.fmu_adresses, tcur, step_size/2.0, False, calling_setting)
1118               
1119                #Exchange and set new inputs
1120                self.set_current_step_size(step_size/2.0)
1121                self.set_input(tcur + step_size/2.0)
1122                ycur, ydcur, ucur = self.exchange_connection_data()
1123                self.set_last_y(ycur)
1124                self.set_last_yd(ydcur)
1125               
1126                #Take another half step
1127                perform_do_step(self.models, self.elapsed_time, self.fmu_adresses, tcur+step_size/2.0, step_size/2.0, False, calling_setting)
1128               
1129                self.exchange_connection_data()
1130                self.set_last_y(ycur)
1131                self.set_last_yd(ydcur)
1132                y_half = self.y_prev.copy()#self.correct_output(self.get_connection_outputs(), y_old)
1133               
1134                step_size_old = step_size
1135               
1136                error = self.estimate_error(y_half, y_full)
1137                step_size = self.adapt_stepsize(step_size, error)
1138                if opts["maxh"] > 0: #Adjust for the maximum step-size (if set)
1139                    step_size = min(step_size, opts["maxh"])
1140               
1141                if opts["logging"]:
1142                    self._error_data["time"].append(tcur)
1143                    self._error_data["error"].append(error)
1144                    self._error_data["step-size"].append(step_size)
1145               
1146                if error < 1.1: #Step accepted
1147                    #Update the time
1148                    tcur += step_size_old
1149                    #Store data
1150                    time_start = timer()
1151                    #store_communication_point(self.models_dict)
1152                    self.report_solution(tcur)
1153                    self.elapsed_time["result_handling"] += timer() - time_start
1154               
1155                    self.statistics["nsteps"] += 1
1156                    y_old = y_half.copy()
1157                   
1158                    if tcur+step_size > final_time: #Make sure that we don't step over the final time
1159                        step_size = final_time - tcur
1160                else:
1161                    #Restore FMU states
1162                    set_fmu_states(states)
1163                    self.y_prev = y_old
1164                    self.set_last_y(y_old)
1165                    self.statistics["nreject"] += 1
1166                   
1167                    if opts["logging"]:
1168                        self._error_data["rejected"].append((tcur, error))
1169           
1170            if states:
1171                free_fmu_states(states)
1172        else:
1173            self.set_current_step_size(step_size)
1174            for tcur in np.arange(start_time, final_time, step_size):
1175                if tcur + step_size > final_time:
1176                    step_size = final_time - tcur
1177                    self.set_current_step_size(step_size)
1178                   
1179                perform_do_step(self.models, self.elapsed_time, self.fmu_adresses, tcur, step_size, True, calling_setting)
1180               
1181                if self.opts["store_step_before_update"]:
1182                    time_start = timer()
1183                    #store_communication_point(self.models_dict)
1184                    self.report_solution(tcur)
1185                    self.elapsed_time["result_handling"] += timer() - time_start
1186                   
1187                #Set external input
1188                self.set_input(tcur + step_size)
1189                ycur, ydcur, ucur = self.exchange_connection_data()
1190                self.set_last_y(ycur)
1191                self.set_last_yd(ydcur)
1192               
1193                time_start = timer()
1194                #store_communication_point(self.models_dict)
1195                self.report_solution(tcur)
1196                self.elapsed_time["result_handling"] += timer() - time_start
1197               
1198                self.statistics["nsteps"] += 1
1199               
1200                #Logging
1201                if opts["logging"]:
1202                    D = self.compute_global_D()
1203                    DL = D.dot(self.L)
1204                    import numpy.linalg
1205                    import scipy.linalg as slin
1206                    print("At time: %E , rho(DL)=%s"%(tcur + step_size, str(numpy.linalg.eig(DL.todense())[0])))
1207                    C = self.compute_global_C()
1208                    if C is not None:
1209                        C = C.todense()
1210                        I = np.eye(*DL.shape)
1211                        LIDLC = self.L.dot(lin.solve(I-DL,C))
1212                        print("           , rho(L(I-DL)^(-1)C)=%s"%(str(numpy.linalg.eig(LIDLC)[0])))
1213                    A = self.compute_global_A()
1214                    B = self.compute_global_B()
1215                    if C is not None and A is not None and B is not None:
1216                        A = A.todense(); B = B.todense()
1217                        eAH = slin.expm(A*step_size)
1218                        K1  = lin.solve(I-DL,C)
1219                        K2  = lin.solve(A,(eAH-np.eye(*eAH.shape)).dot(B.dot(self.L.todense())))
1220                        R1  = np.hstack((eAH, K1))
1221                        R2  = np.hstack((K2.dot(eAH), K2.dot(K1)))
1222                        G   = np.vstack((R1,R2))
1223                        G1  = K2.dot(K1)
1224                        print("           , rho(G)=%s"%(str(numpy.linalg.eig(G1)[0])))
1225                   
1226   
1227    def specify_external_input(self, input):
1228        input_names = input[0]
1229        if isinstance(input_names,tuple):
1230            input_names = [input_names]
1231           
1232        if hasattr(input[1],"__call__"):
1233            input_traj=(input_names,
1234                    TrajectoryUserFunction(input[1]))
1235        else:
1236            input_traj=(input_names,
1237                    TrajectoryLinearInterpolation(input[1][:,0],
1238                                                  input[1][:,1:]))
1239        self.input_traj = input_traj
1240       
1241    cpdef set_input(self, double time):
1242        if self.input_traj is not None:
1243            u = self.input_traj[1].eval(np.array([time]))[0,:]
1244           
1245            for i,m in enumerate(self.input_traj[0]):
1246                m[0].set(m[1],u[i])
1247   
1248    cdef double estimate_error(self, np.ndarray y_half, np.ndarray y_full, int order=0):
1249        cdef np.ndarray err = np.abs((y_half - y_full)/(1-2**(order+1)))
1250       
1251        return np.sqrt(1.0/len(y_half)*sum(err/(self.atol+self.rtol*abs(y_half))))
1252   
1253    cdef double adapt_stepsize(self, double step_size, double error, int order=0):
1254        """
1255        Adjust the stepsize depending on the error.
1256        """
1257        #Get the extrapolation order
1258        cdef double one_over_p = 1.0/(order+2.0)
1259        cdef double alpha = 0.9
1260        cdef double fac1 = 6.0
1261        cdef double fac2 = 0.2
1262       
1263        if error == 0.0:
1264            return step_size*fac1
1265        else:
1266            return step_size*min(fac1,max(fac2,alpha*(1.0/error)**one_over_p))
1267           
1268    cdef reset_statistics(self):
1269        for key in self.elapsed_time: 
1270            self.elapsed_time[key] = 0.0
1271        for key in self.elapsed_time_init: 
1272            self.elapsed_time_init[key] = 0.0
1273        for key in self.statistics.keys():
1274            self.statistics[key] = 0
1275   
1276    def simulate_options(self):
1277        opts = MasterAlgOptions(self)
1278
1279        return opts
1280       
1281    def simulate_profile(self, double start_time=0.0, double final_time=1.0, input=None, options={}):
1282        import pstats, cProfile
1283        res = None
1284        cProfile.runctx("res = self.simulate(start_time, final_time, input, options)", globals(), locals(), "Profile.prof")
1285        s = pstats.Stats("Profile.prof")
1286        s.strip_dirs().sort_stats("cumulative").print_stats(30)
1287       
1288        return self.get_results(options)
1289   
1290    def simulate(self, double start_time=0.0, double final_time=1.0, input=None, options={}):
1291        """
1292        Simulates the system.
1293        """
1294        self.reset_statistics()
1295       
1296        if isinstance(options, dict) and not \
1297            isinstance(options, MasterAlgOptions):
1298            # user has passed dict with options or empty dict = default
1299            options = MasterAlgOptions(self, options)
1300        elif isinstance(options, MasterAlgOptions):
1301            pass
1302        else:
1303            raise InvalidAlgorithmOptionException(options)
1304       
1305        #if options == "Default":
1306        #   options = self.simulate_options()
1307        if options["linear_correction"]:
1308            if self.support_directional_derivatives != 1 and not options["experimental_finite_difference_D"]:
1309                warnings.warn("Linear correction only supported if directional derivatives are available.")
1310                self.linear_correction = 0
1311            else:
1312                self.linear_correction = 1
1313        else:
1314            self.linear_correction = 0
1315        if options["error_controlled"]:
1316            if self.support_storing_fmu_states != 1:
1317                warnings.warn("Error controlled simulation only supported if storing FMU states are available.")
1318                self.error_controlled = 0
1319            else:
1320                self.error_controlled = 1
1321                self.atol = options["atol"]
1322                self.rtol = options["rtol"]
1323        else:
1324            self.error_controlled = 0
1325        if options["extrapolation_order"] > 0:
1326            if self.support_interpolate_inputs != 1:
1327                warnings.warn("Extrapolation of inputs only supported if the individual FMUs support interpolation of inputs.")
1328                options["extrapolation_order"] = 0
1329       
1330        if options["num_threads"] and options["execution"] == "parallel":
1331            pass
1332            IF WITH_OPENMP: 
1333                openmp.omp_set_num_threads(options["num_threads"])
1334        if options["step_size"] <= 0.0:
1335            raise fmi.FMUException("The step-size must be greater than zero.")
1336        if options["maxh"] < 0.0:
1337            raise fmi.FMUException("The maximum step-size must be greater than zero (or equal to zero, i.e inactive).")
1338       
1339        self.opts = options #Store the options
1340           
1341        if input is not None:
1342            self.specify_external_input(input)
1343       
1344        #Initialize the models
1345        if options["initialize"]:
1346            time_start = timer()
1347            self.initialize(start_time, final_time, options)
1348            time_stop  = timer()
1349            print('Elapsed initialization time: ' + str(time_stop-time_start) + ' seconds.')
1350       
1351        #Store the inputs
1352        self.set_last_us(self.L.dot(self.get_connection_outputs()))
1353       
1354        #Copy FMU address (used when evaluating in parallel)
1355        self.copy_fmu_addresses()
1356       
1357        self.initialize_result_objects(options)
1358        store_communication_point(self.models_dict)
1359        #self.report_solution(tcur)
1360       
1361        #Start of simulation, start the clock
1362        #time_start = time.clock()
1363        time_start = timer()
1364        self._time_integration_start = time_start
1365       
1366        self.jacobi_algorithm(start_time,final_time, options)
1367
1368        #End of simulation, stop the clock
1369        #time_stop = time.clock()
1370        time_stop = timer()
1371       
1372        self.print_statistics(options)
1373        print('')
1374        print('Simulation interval      : ' + str(start_time) + ' - ' + str(final_time) + ' seconds.')
1375        print('Elapsed simulation time  : ' + str(time_stop-time_start) + ' seconds.')
1376        for model in self.models:
1377            print(' %f seconds spent in %s.'%(self.elapsed_time[model],model.get_name()))
1378        print(' %f seconds spent saving simulation result.'%(self.elapsed_time["result_handling"]))
1379       
1380        #Write the results to file and return
1381        return self.get_results(options)
1382   
1383    def get_results(self, opts):
1384        """
1385        For each model, write the result and load the result file
1386        """
1387        finalize_result_objects(self.models_dict)
1388       
1389        res = {}
1390        #Load data
1391        for i,model in enumerate(self.models):
1392            if opts["result_handling"] == "file":
1393                dym_textual = ResultDymolaTextual(model.get_identifier()+'_'+str(i)+'_result.txt')
1394                res[i] = AssimuloSimResult(model, model.get_identifier()+'_'+str(i)+'_result.txt', None, dym_textual, None)
1395                res[model] = res[i]
1396            elif opts["result_handling"] == "none":
1397                res[model] = None
1398            else:
1399                raise fmi.FMUException("Currently only writing result to file or none is supported.")
1400       
1401        return res
1402   
1403    def print_statistics(self, opts):
1404        print('Master Algorithm options:')
1405        print(' Algorithm             : Jacobi ' + ("(variable-step)" if self.error_controlled else "(fixed-step)"))
1406        print('  Execution            : ' + ("Parallel" if self.opts["execution"] == "parallel" else "Serial"))
1407        print(' Extrapolation Order   : ' + str(opts["extrapolation_order"]) + ("(with smoothing)" if opts["smooth_coupling"] and opts["extrapolation_order"] > 0  else ""))
1408        if self.error_controlled:
1409            print(' Tolerances (relative) : ' + str(opts["rtol"]))
1410            print(' Tolerances (absolute) : ' + str(opts["atol"]))
1411        else:
1412            print(' Step-size             : ' + str(opts["step_size"]))
1413        print(' Algebraic loop        : ' + ("True" if self.algebraic_loops else "False"))
1414        print('  Linear Correction    : ' + ("True" if self.linear_correction else "False"))
1415        print('')
1416        print('Statistics: ')
1417        print(' Number of global steps        : %d'%self.statistics["nsteps"])
1418        if self.error_controlled:
1419            print(' Number of rejected steps      : %d'%self.statistics["nreject"])
1420   
1421    def _get_support_directional_derivatives(self):
1422        if self._support_directional_derivatives == -1:
1423            self._support_directional_derivatives = 1
1424            for model in self.models_dict.keys():
1425                if model.get_capability_flags()["providesDirectionalDerivatives"] is False:
1426                    self._support_directional_derivatives = 0
1427                    break 
1428        return self._support_directional_derivatives
1429
1430    support_directional_derivatives = property(_get_support_directional_derivatives)
1431   
1432    def _get_support_storing_fmu_states(self):
1433        if self._support_storing_fmu_states == -1:
1434            self._support_storing_fmu_states = 1
1435            for model in self.models_dict.keys():
1436                if model.get_capability_flags()["canGetAndSetFMUstate"] is False:
1437                    self._support_storing_fmu_states = 0
1438                    break
1439        return self._support_storing_fmu_states
1440   
1441    support_storing_fmu_states = property(_get_support_storing_fmu_states)
1442   
1443    def _get_support_interpolate_inputs(self):
1444        if self._support_interpolate_inputs == -1:
1445            self._support_interpolate_inputs = 1
1446            for model in self.models_dict.keys():
1447                if model.get_capability_flags()["canInterpolateInputs"] is False:
1448                    self._support_interpolate_inputs = 0
1449                    break
1450        return self._support_interpolate_inputs
1451
1452    support_interpolate_inputs = property(_get_support_interpolate_inputs)
1453   
1454    def _get_max_output_derivative_order(self):
1455        if self._max_output_derivative_order == -1:
1456            self._max_output_derivative_order = 10
1457            for model in self.models_dict.keys():
1458                self._max_output_derivative_order = min(self._max_output_derivative_order, model.get_capability_flags()["maxOutputDerivativeOrder"])
1459        return self._max_output_derivative_order
1460   
1461    max_output_derivative_order = property(_get_max_output_derivative_order)
1462   
1463    def reset(self):
1464        """
1465        Reset the coupled models.
1466        """
1467        reset_models(self.models)
1468   
1469    def visualize_connections(self, vectorized=True, delete=False):
1470        import subprocess
1471        import tempfile
1472        import os
1473       
1474        def strip_var(var): return var.replace(".","_").replace("[","_").replace("]","_").replace(",","_")
1475        def vec_vars(list_vars):
1476            if not isinstance(list_vars, list): list_vars = [list_vars]
1477            res_vars = OrderedDict()
1478            for var in list_vars:
1479                if fnmatch.fnmatch(var, "*[[]?[]]"):
1480                    res_vars[var[:-3]] = 1
1481                elif fnmatch.fnmatch(var, "*[[]?,?[]]"):
1482                    res_vars[var[:-5]] = 1
1483                else:
1484                    res_vars[var] = 1
1485            return res_vars.keys()
1486       
1487        tmp = next(tempfile._get_candidate_names())
1488       
1489        with open("%s.dot"%tmp, 'w') as dfile:
1490            dfile.write("""digraph "callgraph" {
1491    nodesep=0.5;ranksep=1;
1492    node [fontname=Arial, height=0, shape=box, style=filled, width=0];
1493    edge [fontname=Arial];\n""")
1494            for model in self.models_dict:
1495                #dfile.write('   subgraph cluster_%s { fontsize="12"; label = "%s"; \n'%(model.get_name().replace(".","_"), model.get_name()))
1496                dfile.write('   subgraph cluster_%s { fontsize="12"; label = "%s"; \n'%(str(id(model)), model.get_name()))
1497               
1498                if vectorized:
1499                    input_vars = vec_vars(self.models_dict[model]["local_input"])
1500                else:
1501                    input_vars = self.models_dict[model]["local_input"]
1502               
1503                tmp_str = "<input> Inputs|"
1504                for input_var in input_vars:
1505                    tmp_str += "<%s> %s|"%(strip_var(input_var), input_var) 
1506                dfile.write('       input_%s[label="{%s}" shape=Mrecord];\n'%(str(id(model)), tmp_str[:-1]))
1507               
1508                if vectorized:
1509                    output_vars = vec_vars(self.models_dict[model]["local_output"])
1510                else:
1511                    output_vars = self.models_dict[model]["local_output"]
1512               
1513                tmp_str = "<output> Outputs|"
1514                for output_var in output_vars:
1515                    tmp_str += "<%s> %s|"%(strip_var(output_var), output_var)
1516                dfile.write('       output_%s[label="{%s}" shape=Mrecord]; }\n'%(str(id(model)), tmp_str[:-1]))
1517                #dfile.write("   }\n")
1518               
1519                for direct_dependence in self.models_dict[model]["direct_dependence"]:
1520                    if vectorized:
1521                        dfile.write("   input_%s:%s:e -> output_%s:%s:w[color=red]; \n"%(str(id(model)),strip_var(vec_vars(direct_dependence[0])[0]),str(id(model)),strip_var(vec_vars(direct_dependence[1])[0])))
1522                    else:
1523                        dfile.write("   input_%s:%s:e -> output_%s:%s:w[color=red]; \n"%(str(id(model)),strip_var(direct_dependence[0]),str(id(model)),strip_var(direct_dependence[1])))
1524
1525            for connection in self.connections:
1526                if vectorized:
1527                    dfile.write("   output_%s:%s:e -> input_%s:%s:w; \n"%(str(id(connection[0])),strip_var(vec_vars(connection[1])[0]),str(id(connection[2])),strip_var(vec_vars(connection[3])[0])))
1528                else:
1529                    dfile.write("   output_%s:%s:e -> input_%s:%s:w; \n"%(str(id(connection[0])),strip_var(connection[1]),str(id(connection[2])),strip_var(connection[3])))
1530            dfile.write("}\n")
1531       
1532        with open("%s.dot"%tmp, 'r') as dfile:
1533            lines = (line.rstrip() for line in dfile)
1534            unique_lines = OrderedDict.fromkeys((line for line in lines if line))
1535        with open("%s.dot"%tmp, 'w') as dfile:
1536            for line in unique_lines:
1537                dfile.write(line+"\n")
1538       
1539        subprocess.call(["dot -Tps %s.dot -o %s.pdf"%(tmp,tmp)], shell=True)
1540        #subprocess.Popen([os.path.join(os.getcwd(),"%s.pdf"%tmp)],shell=True)
1541        os.system("/usr/bin/xdg-open %s.pdf"%tmp)
1542       
1543    def define_graph(self):
1544        edges = []
1545        graph_info = OrderedDict()
1546       
1547        def add_info(primary_key, secondary_key, info):
1548            try:
1549                if isinstance(info, list):
1550                    graph_info[primary_key][secondary_key].extend(info)
1551                else:
1552                    graph_info[primary_key][secondary_key] = info
1553            except KeyError:
1554                graph_info[primary_key] = {secondary_key: info}
1555       
1556        #Define the edges:
1557        for model in self.models:
1558            for direct_dependence in self.models_dict[model]["direct_dependence"]:
1559                add_info(str(id(model))+"_"+direct_dependence[0], "type", 0) #Input
1560                add_info(str(id(model))+"_"+direct_dependence[1], "type", 1) #Output
1561                add_info(str(id(model))+"_"+direct_dependence[0], "model", str(id(model)))
1562                add_info(str(id(model))+"_"+direct_dependence[1], "model", str(id(model)))
1563                add_info(str(id(model)), "outputs", [str(id(model))+"_"+direct_dependence[1]])
1564               
1565                edges.append((str(id(model))+"_"+direct_dependence[0],str(id(model))+"_"+direct_dependence[1]))
1566        for connection in self.connections:
1567            add_info(str(id(connection[0]))+"_"+connection[1], "type", 1) #Output
1568            add_info(str(id(connection[2]))+"_"+connection[3], "type", 0) #Input
1569            add_info(str(id(connection[0]))+"_"+connection[1], "model", str(id(connection[0]))) #Output
1570            add_info(str(id(connection[2]))+"_"+connection[3], "model", str(id(connection[2]))) #Input
1571            add_info(str(id(connection[0])), "outputs", [str(id(connection[0]))+"_"+connection[1]])
1572           
1573            edges.append((str(id(connection[0]))+"_"+connection[1],str(id(connection[2]))+"_"+connection[3]))
1574       
1575        return Graph(edges, graph_info)
1576   
1577    def compute_evaluation_order(self, init_type="greedy", order = None):
1578       
1579        if order is None:
1580            graph = self.define_graph()
1581           
1582            if init_type == "grouping":
1583                order = graph.grouped_order(graph.strongly_connected_components())[::-1]
1584            elif init_type == "simple":
1585                order = graph.strongly_connected_components()[::-1]
1586            elif init_type == "greedy":
1587                order = graph.compute_evaluation_order()[::-1]
1588            else:
1589                raise Exception("Unknown initialization type. Use either 'greedy', 'simple', 'grouping'")
1590       
1591        blocks = []
1592        for block in order:
1593            blocks.append({})
1594            blocks[-1]["inputs"] = {}
1595            blocks[-1]["outputs"] = {}
1596            blocks[-1]["has_outputs"] = False
1597            for variable_compound in block:
1598                model_id, var = variable_compound.split("_",1)
1599                model = self.models_id_mapping[model_id]
1600                if var in self.models_dict[model]["local_input"]:
1601                    try:
1602                        blocks[-1]["inputs"][model].append(var)
1603                    except KeyError:
1604                        blocks[-1]["inputs"][model] = [var]
1605                elif var in self.models_dict[model]["local_output"]:
1606                    try:
1607                        blocks[-1]["outputs"][model].append(var)
1608                    except KeyError:
1609                        blocks[-1]["outputs"][model] = [var]
1610                        blocks[-1]["has_outputs"] = True
1611                else:
1612                    raise fmi.FMUException("Something went wrong while creating the blocks.")
1613       
1614        for block in blocks:
1615            block["global_inputs_mask"] = np.array([False]*self._len_inputs)
1616            block["global_outputs_mask"] = np.array([False]*self._len_outputs)
1617            block["inputs_mask"] = {}
1618            block["outputs_mask"] = {}
1619            for model in block["inputs"].keys():
1620                input_vref = np.array(self.models_dict[model]["local_input"])
1621                mask = np.array([False]*len(input_vref))
1622                for i,x in enumerate(block["inputs"][model]):
1623                    pos = np.where(input_vref == x)[0][0]
1624                    mask[pos] = True
1625                    block["global_inputs_mask"][self.models_dict[model]["global_index_inputs"]+pos] = True
1626                block["inputs_mask"][model] = mask
1627                #if len(block["inputs"].keys()) > 1:
1628                #    print "Creation: ", model, block["inputs_mask"]
1629                #print block["inputs"].keys(), model, block["inputs_mask"]
1630            for model in block["outputs"].keys():
1631                output_vref = np.array(self.models_dict[model]["local_output"])
1632                mask = np.array([False]*len(output_vref))
1633                for i,x in enumerate(block["outputs"][model]):
1634                    pos = np.where(output_vref == x)[0][0]
1635                    mask[pos] = True
1636                    block["global_outputs_mask"][self.models_dict[model]["global_index_outputs"]+pos] = True
1637                block["outputs_mask"][model] = mask
1638       
1639        #Compress blocks
1640        compressed_blocks = []
1641        """
1642        compressed_blocks = [blocks[0]]
1643        last_has_inputs = True if len(blocks[0]["inputs"]) > 0 else False
1644        last_has_outputs = True if len(blocks[0]["outputs"]) > 0 else False
1645        for block in blocks[1:]:
1646            has_inputs = True if len(block["inputs"]) > 0 else False
1647            has_outputs = True if len(block["outputs"]) > 0 else False
1648            if has_inputs and has_outputs or last_has_inputs and last_has_outputs:
1649                compressed_blocks.append(block)
1650                last_has_inputs = has_inputs
1651                last_has_outputs = has_outputs
1652            elif last_has_inputs and not last_has_outputs and has_inputs and not has_outputs:
1653                for key in block["inputs"].keys():
1654                    try:
1655                        compressed_blocks[-1]["inputs"][key].append(block["inputs"][key])
1656                    except KeyError:
1657                        compressed_blocks[-1]["inputs"][key] = block["inputs"][key]
1658            elif not last_has_inputs and last_has_outputs and not has_inputs and has_outputs:
1659                for key in block["outputs"].keys():
1660                    try:
1661                        compressed_blocks[-1]["outputs"][key].append(block["outputs"][key])
1662                    except KeyError:
1663                        compressed_blocks[-1]["outputs"][key] = block["outputs"][key]
1664            else:
1665                compressed_blocks.append(block)
1666                last_has_inputs = has_inputs
1667                last_has_outputs = has_outputs
1668        """
1669        return order, blocks,compressed_blocks
1670       
1671       
Note: See TracBrowser for help on using the repository browser.