diff --git a/src/oemof/solph/__init__.py b/src/oemof/solph/__init__.py
index 8cf2b6cb1..0532ecd5c 100644
--- a/src/oemof/solph/__init__.py
+++ b/src/oemof/solph/__init__.py
@@ -10,6 +10,7 @@
 from .components import OffsetTransformer  # noqa: F401
 from .groupings import GROUPINGS  # noqa: F401
 from .models import Model  # noqa: F401
+from .models import MultiObjectiveModel  # noqa: F401
 from .network import Bus  # noqa: F401
 from .network import EnergySystem  # noqa: F401
 from .network import Flow  # noqa: F401
@@ -17,6 +18,7 @@
 from .network import Source  # noqa: F401
 from .network import Transformer  # noqa: F401
 from .options import Investment  # noqa: F401
+from .options import MultiObjective  # noqa: F401
 from .options import NonConvex  # noqa: F401
 from .plumbing import sequence  # noqa: F401
 from .processing import parameter_as_dict  # noqa: F401
diff --git a/src/oemof/solph/blocks/__init__.py b/src/oemof/solph/blocks/__init__.py
index bf7dc70b9..56699a17a 100644
--- a/src/oemof/solph/blocks/__init__.py
+++ b/src/oemof/solph/blocks/__init__.py
@@ -7,5 +7,6 @@
 from .bus import Bus  # noqa: F401
 from .flow import Flow  # noqa: F401
 from .investment_flow import InvestmentFlow  # noqa: F401
+from .multiobjective_flow import MultiObjectiveFlow  # noqa: F401
 from .non_convex_flow import NonConvexFlow  # noqa: F401
 from .transformer import Transformer  # noqa: F401
diff --git a/src/oemof/solph/blocks/multiobjective_flow.py b/src/oemof/solph/blocks/multiobjective_flow.py
new file mode 100644
index 000000000..22478dfb4
--- /dev/null
+++ b/src/oemof/solph/blocks/multiobjective_flow.py
@@ -0,0 +1,92 @@
+# -*- coding: utf-8 -*-
+
+"""Creating sets, variables, constraints and parts of the objective function
+for MultiObjectiveFlow objects.
+"""
+
+from collections import defaultdict
+
+from pyomo.core import Set
+from pyomo.core.base.block import SimpleBlock
+
+
+class MultiObjectiveFlow(SimpleBlock):
+    r""" Block for all flows with :attr:`multiobjective` being not None.
+
+    Adds no additional variables or constraints. Used solely to add index set
+    and objective expressions.
+
+    **The following sets are created:** (-> see basic sets at :class:`.Model` )
+
+    MULTIOBJECTIVEFLOWS
+        A set of flows with the attribute :attr:`multiobjective` of type
+        :class:`.options.MultiObjective`.
+
+    **The following parts of the objective function are created:**
+
+    If :attr:`variable_costs` are set by the user:
+      .. math::
+          \sum_{(i,o)} \sum_t flow(i, o, t) \cdot variable\_costs(e, i, o, t)
+
+          \forall e \in \{set\ of\ objective\ functions\}
+
+    Specific objectives are created in the same manner as for :class:`Flow`.
+
+    """
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+
+    def _create(self, group=None):
+        r"""Creates sets, variables and constraints for all standard flows.
+
+        Parameters
+        ----------
+        group : list
+            List containing tuples containing flow (f) objects and the
+            associated source (s) and target (t)
+            of flow e.g. groups=[(s1, t1, f1), (s2, t2, f2),..]
+        """
+        if group is None:
+            return None
+
+        # ########################## SETS #################################
+        self.MULTIOBJECTIVEFLOWS = Set(
+            initialize=[(g[0], g[1]) for g in group])
+
+    def _objective_expression(self):
+        r"""Objective expression for all multiobjective flows with fixed costs
+        and variable costs.
+        """
+        if not hasattr(self, "MULTIOBJECTIVEFLOWS"):
+            return 0
+        m = self.parent_block()
+
+        mo_costs = defaultdict(lambda: 0, {'_standard': 0})
+
+        for i, o in self.MULTIOBJECTIVEFLOWS:
+            for mo_key, mo_param in m.flows[i, o].multiobjective.items():
+                variable_costs = 0
+                gradient_costs = 0
+                if mo_param.variable_costs[0] is not None:
+                    for t in m.TIMESTEPS:
+                        variable_costs += (
+                            m.flow[i, o, t]
+                            * m.objective_weighting[t]
+                            * mo_param.variable_costs[t])
+
+                if mo_param.positive_gradient["ub"][0] is not None:
+                    for t in m.TIMESTEPS:
+                        gradient_costs += (
+                            self.positive_gradient[i, o, t]
+                            * mo_param.positive_gradient["costs"])
+
+                if mo_param.negative_gradient["ub"][0] is not None:
+                    for t in m.TIMESTEPS:
+                        gradient_costs += (
+                            self.negative_gradient[i, o, t]
+                            * mo_param.negative_gradient["costs"])
+
+                mo_costs[mo_key] += variable_costs + gradient_costs
+
+        return mo_costs
diff --git a/src/oemof/solph/groupings.py b/src/oemof/solph/groupings.py
index bfd71138d..d11624864 100644
--- a/src/oemof/solph/groupings.py
+++ b/src/oemof/solph/groupings.py
@@ -87,9 +87,24 @@ def _nonconvex_grouping(stf):
 )
 
 
+def _multiobjective_grouping(stf):
+    if hasattr(stf[2], "multiobjective"):
+        if stf[2].multiobjective is not None:
+            return True
+    else:
+        return False
+
+
+multiobjective_flow_grouping = groupings.FlowsWithNodes(
+    constant_key=blocks.MultiObjectiveFlow,
+    filter=_multiobjective_grouping
+)
+
+
 GROUPINGS = [
     constraint_grouping,
     investment_flow_grouping,
     standard_flow_grouping,
     nonconvex_flow_grouping,
+    multiobjective_flow_grouping,
 ]
diff --git a/src/oemof/solph/models.py b/src/oemof/solph/models.py
index 04bf0e0f8..dc9c5292c 100644
--- a/src/oemof/solph/models.py
+++ b/src/oemof/solph/models.py
@@ -13,6 +13,10 @@
 """
 import logging
 import warnings
+from collections import defaultdict
+
+import numpy as np
+from pandas import DataFrame
 
 from pyomo import environ as po
 from pyomo.core.plugins.transform.relax_integrality import RelaxIntegrality
@@ -362,3 +366,323 @@ def _add_parent_block_variables(self):
                 if (o, i) in self.UNIDIRECTIONAL_FLOWS:
                     for t in self.TIMESTEPS:
                         self.flow[o, i, t].setlb(0)
+
+
+class MultiObjectiveModel(Model):
+    """An  energy system model for operational and investment
+    optimization.
+
+    Parameters
+    ----------
+    energysystem : EnergySystem object
+        Object that holds the nodes of an oemof energy system graph
+    constraint_groups : list
+        Solph looks for these groups in the given energy system and uses them
+        to create the constraints of the optimization problem.
+        Defaults to `Model.CONSTRAINTS`
+
+    **The following basic sets are created**:
+
+    NODES :
+        A set with all nodes of the given energy system.
+
+    TIMESTEPS :
+        A set with all timesteps of the given time horizon.
+
+    FLOWS :
+        A 2 dimensional set with all flows. Index: `(source, target)`
+
+    **The following basic variables are created**:
+
+    flow
+        Flow from source to target indexed by FLOWS, TIMESTEPS.
+        Note: Bounds of this variable are set depending on attributes of
+        the corresponding flow object.
+
+    """
+
+    CONSTRAINT_GROUPS = [
+        blocks.Bus,
+        blocks.Transformer,
+        blocks.InvestmentFlow,
+        blocks.Flow,
+        blocks.NonConvexFlow,
+        blocks.MultiObjectiveFlow
+    ]
+
+    def __init__(self, energysystem, **kwargs):
+        super().__init__(energysystem, **kwargs)
+
+    def _add_objective(self, sense=po.minimize, update=False):
+        """ Method to sum up all objective expressions from the child blocks
+        that have been created. This method looks for `_objective_expression`
+        attribute in the block definition and will call this method to add
+        their return value to the objective function or to the respective
+        objective function if multiple objective function keys are given.
+        """
+
+        if update:
+            self.del_component('objective')
+
+        # create dict for all distinct objective expressions
+        self.objective_functions = defaultdict(lambda: 0, {'_standard': 0})
+
+        # set sign based on optimisation sense
+        if sense == po.minimize:
+            sign = 1
+        else:
+            sign = -1
+
+        for block in self.component_data_objects():
+            if hasattr(block, '_objective_expression'):
+                expr = block._objective_expression()
+                if isinstance(expr, defaultdict):
+                    for obj_key, obj_val in expr.items():
+                        self.objective_functions[obj_key] += (
+                            sign * obj_val)
+                else:
+                    self.objective_functions['_standard'] += (
+                        sign * expr)
+
+    def solve(self, solver='cbc', solver_io='lp', **kwargs):
+        r""" Takes care of communication with solver to solve the model.
+        Differentiates between single objective optimization or multi
+        objective optimization. Defaults to single objective optimization.
+
+        Parameters
+        ----------
+        solver : string
+            solver to be used e.g. "glpk","gurobi","cplex"
+        solver_io : string
+            pyomo solver interface file format: "lp","python","nl", etc.
+        \**kwargs : keyword arguments
+            Possible keys can be set see below:
+
+        Other Parameters
+        ----------------
+        solve_kwargs : dict
+            Other arguments for the pyomo.opt.SolverFactory.solve() method
+            Example : {"tee":True}
+        cmdline_options : dict
+            Dictionary with command line options for solver e.g.
+            {"mipgap":"0.01"} results in "--mipgap 0.01"
+            {"interior":" "} results in "--interior"
+            Gurobi solver takes numeric parameter values such as
+            {"method": 2}
+        optimization_type: str, default 'singular'
+            Sets type of optimization, currently either 'singular' for single
+            objective or 'weighted' for weighted sum of several objectives.
+        objective (with 'singular'): str, default '_standard'
+            Name of singular objective function to use. Defaults to internal
+            standard name.
+        objective_weights (with 'weighted'): dict, default {'_standard': 1}
+            Dictionary with names of objectives as keys and numeric weights
+            as values.
+        """
+        solve_kwargs = kwargs.get('solve_kwargs', {})
+        solver_cmdline_options = kwargs.get("cmdline_options", {})
+
+        opt = SolverFactory(solver, solver_io=solver_io)
+        # set command line options
+        options = opt.options
+        for k in solver_cmdline_options:
+            options[k] = solver_cmdline_options[k]
+
+        # get type of optimisation
+        optimization_type = kwargs.get('optimization_type', 'singular')
+
+        # delete objective if it exists
+        self.del_component('objective')
+
+        if optimization_type == 'singular':
+            # get function to use
+            objective = kwargs.get('objective', '_standard')
+
+            if not isinstance(objective, str):
+                raise TypeError('Objective is not of type "string"')
+            if objective not in self.objective_functions.keys():
+                raise ValueError(
+                    'No cost for objective "{0}"'.format(objective))
+
+            # log objective
+            logging.info('Objective set to {0}.'.format(objective))
+
+            # set chosen objective function
+            self.objective = po.Objective(
+                sense=po.minimize,
+                expr=self.objective_functions.get(objective, 0.0))
+
+            solver_results = opt.solve(self, **solve_kwargs)
+
+            status = solver_results["Solver"][0]["Status"]
+            termination_condition = (
+                solver_results["Solver"][0]["Termination condition"])
+
+            if status == "ok" and termination_condition == "optimal":
+                logging.info("Optimization successful...")
+            else:
+                msg = ("Optimization ended with status {0} and termination "
+                       "condition {1}")
+                warnings.warn(msg.format(status, termination_condition),
+                              UserWarning)
+            self.es.results = solver_results
+            self.solver_results = solver_results
+
+            return solver_results
+
+        elif optimization_type == 'weighted':
+
+            # get function names and weights to use
+            obj_weights = kwargs.get(
+                'objective_weights',
+                {'_standard': 1})
+
+            # check for correct type
+            if not isinstance(obj_weights, dict):
+                raise TypeError("Objective weights must be of type 'dict'.")
+
+            # check existance of objectives
+            if len(obj_weights) == 0:
+                raise ValueError("Objective weights must not be empty.")
+
+            # initialise weighted sum of objectives
+            expr = 0
+
+            # check and set objectives and weights
+            for obj_name, obj_weight in obj_weights.items():
+                if not isinstance(obj_name, str):
+                    raise TypeError(
+                        'Objective "{0}" is not of type "string"'.format(
+                            obj_name))
+                if obj_name not in self.objective_functions.keys():
+                    raise ValueError('No cost for objective "{0}"'.format(
+                        obj_name))
+
+                # add to summed objective
+                expr += (obj_weight
+                         * self.objective_functions.get(obj_name))
+
+            # log objective and weights
+            logging.info('Objectives set to: {0} with weights: {1}'.format(
+                ', '.join(obj_weights.keys()), str(obj_weights)))
+
+            # create objective
+            self.objective = po.Objective(sense=po.minimize, expr=expr)
+
+            # solve
+            solver_results = opt.solve(self, **solve_kwargs)
+
+            status = solver_results["Solver"][0]["Status"]
+            termination_condition = (
+                solver_results["Solver"][0]["Termination condition"])
+
+            if status == "ok" and termination_condition == "optimal":
+                logging.info("Optimization successful...")
+            else:
+                msg = ("Optimization ended with status {0} and termination "
+                       "condition {1}")
+                warnings.warn(msg.format(status, termination_condition),
+                              UserWarning)
+            self.es.results = solver_results
+            self.solver_results = solver_results
+
+            return solver_results
+        else:
+            raise Exception('Invalid optimization type')
+
+    def pareto(self, solver='cbc', solver_io='lp', **kwargs):
+        solve_kwargs = kwargs.get('solve_kwargs', {})
+        solver_cmdline_options = kwargs.get("cmdline_options", {})
+
+        opt = SolverFactory(solver, solver_io=solver_io)
+        # set command line options
+        options = opt.options
+        for k in solver_cmdline_options:
+            options[k] = solver_cmdline_options[k]
+
+
+        # number of dimensions
+        obj_list = kwargs.get('objectives', list())
+
+        # check if all objectives ar unique
+        objectives = []
+        for obj in obj_list:
+            if obj not in objectives:
+                objectives.append(obj)
+            else:
+                msg = ("Objectives should be given only once."
+                       + " Truncated {{obj_list}} to contain only"
+                       + " unique values")
+                warnings.warn(msg, UserWarning)
+
+        if len(objectives) <= 1:
+            raise ValueError(
+                'List of objectives must contain at'
+                + ' least 2 unique entries!')
+
+        # number of dimensions for grid
+        ndim = len(objectives)
+
+        # number of points per dimension
+        npoints = kwargs.get('npoints', 2)
+
+        # create unscaled grid of weights
+        weights_unscaled = np.meshgrid(
+            *(range(npoints) for i in range(ndim)))
+
+        # reshape (flatten) to 1d arrays
+        weights_unscaled = np.squeeze(np.array(
+            [np.reshape(w, (-1, 1))[1:, :] for w in weights_unscaled]))
+
+        # calculate unique weights scaled to 1
+        weights_scaled = np.unique(
+            weights_unscaled / weights_unscaled.sum(axis=0), axis=1).T
+
+        # create DataFrame to hold results for each point/objective
+        results = DataFrame(
+            index=range(weights_scaled.shape[0]),
+            data=0,
+            columns=objectives)
+
+        # add counter for current weight combination
+        j=0
+
+        # iterate over different weights
+        for weights in weights_scaled:
+            # delete objective if it exists
+            self.del_component('objective')
+
+            # initialise weighted sum of objectives
+            expr = 0
+
+            # iterate over unique objectives
+            for i in range(ndim):
+                # add to summed objective
+                expr += (
+                    weights[i]
+                    * self.objective_functions.get(objectives[i]))
+
+            # create objective
+            self.objective = po.Objective(sense=po.minimize, expr=expr)
+
+            # solve
+            opt.solve(self, **solve_kwargs)
+
+            for i in range(ndim):
+                results.at[j, objectives[i]] = (po.value(
+                    self.objective_functions.get(objectives[i])))
+
+            j += 1
+
+        # find weight combinations belonging to pareto frontier
+        costs = np.array(results)
+        pareto_indices = np.ones(costs.shape[0], dtype = bool)
+        for i, c in enumerate(costs):
+            if pareto_indices[i]:
+                # keep any point with a lower cost
+                pareto_indices[pareto_indices] = (
+                    np.any(costs[pareto_indices]<c, axis=1))
+                # and keep self
+                pareto_indices[i] = True
+
+        return results, weights_scaled, pareto_indices
diff --git a/src/oemof/solph/network/flow.py b/src/oemof/solph/network/flow.py
index 99fef9590..5ee24b0cb 100644
--- a/src/oemof/solph/network/flow.py
+++ b/src/oemof/solph/network/flow.py
@@ -92,17 +92,23 @@ class Flow(on.Edge):
         :class:`Flow <oemof.solph.blocks.Flow>`.
         Note: at the moment this does not work if the investment attribute is
         set .
+    multiobjective: :class:`oemof.solph.options.MultiObjective`
+        Object indicating if multiple objectives should be added for this
+        flow. To be used together with
+        :class:`oemof.solph.models.MultiObjectiveModel` to generate
+        multi objective program.
 
     Notes
     -----
     The following sets, variables, constraints and objective parts are created
-     * :py:class:`~oemof.solph.blocks.flow.Flow`
-     * :py:class:`~oemof.solph.blocks.investment_flow.InvestmentFlow`
-        (additionally if Investment object is present)
-     * :py:class:`~oemof.solph.blocks.non_convex_flow.NonConvexFlow`
-        (If nonconvex  object is present, CAUTION: replaces
-        :py:class:`~oemof.solph.blocks.flow.Flow`
-        class and a MILP will be build)
+     * :py:class:`~oemof.solph.blocks.Flow`
+     * :py:class:`~oemof.solph.blocks.InvestmentFlow` (additionally if
+       Investment object is present)
+     * :py:class:`~oemof.solph.blocks.MultiObjectiveFlow` (additionally if
+       MultiObjective object is present)
+     * :py:class:`~oemof.solph.blocks.NonConvexFlow` (If
+        nonconvex  object is present, CAUTION: replaces
+        :py:class:`~oemof.solph.blocks.Flow` class and a MILP will be build)
 
     Examples
     --------
@@ -137,6 +143,7 @@ def __init__(self, **kwargs):
             "investment",
             "nonconvex",
             "integer",
+            "multiobjective",
         ]
         sequences = ["fix", "variable_costs", "min", "max"]
         dictionaries = ["positive_gradient", "negative_gradient"]
@@ -212,6 +219,20 @@ def __init__(self, **kwargs):
                 + "nonconvex flows!"
             )
 
+        # Checking for attribute combinations with possibly unknown effects
+        if self.multiobjective and sum(self.variable_costs) != 0:
+            msg = (
+                "Multiobjective flows should not be combined with standard"
+                + " variable costs! Costs are automatically attributed to"
+                + " standard objective.")
+            warn(msg, UserWarning)
+        if self.multiobjective and self.investment:
+            msg = (
+                "Multiobjective flows should not be combined with investment"
+                + " flows! Investment is automatically attributed to standard"
+                + " objective.")
+            warn(msg, UserWarning)
+
         # Checking for impossible gradient combinations
         if self.nonconvex:
             if self.nonconvex.positive_gradient["ub"][0] is not None and (
diff --git a/src/oemof/solph/options.py b/src/oemof/solph/options.py
index 3a52dd298..ceade84e4 100644
--- a/src/oemof/solph/options.py
+++ b/src/oemof/solph/options.py
@@ -211,3 +211,77 @@ def max_up_down(self):
             self._calculate_max_up_down()
 
         return self._max_up_down
+
+
+class MultiObjective:
+    """
+    Parameters
+    ----------
+    [key-value-pairs]: MultiObjective.Objective
+        Key value pair consisting of name for partial objective function and
+        instance of Objective nested class with corresponding parameters.
+    """
+
+    def __init__(self, **kwargs):
+
+        # initalise internal dict
+        self.mo = dict()
+
+        for mo_key, mo_param in kwargs.items():
+            if isinstance(mo_param, MultiObjective.Objective):
+                self.mo[mo_key] = mo_param
+
+    def items(self):
+        return self.mo.items()
+
+    class Objective:
+        """
+        Comfort class to faciliate construction of different objectives.
+
+        Comfort class wrapping objective function parameters for singular
+        objective function in multi objective optimization. Uses same
+        parameters for objective declaration as
+        :class:`Flow <oemof.solph.blocks.Flow>`.
+
+        Parameters
+        ----------
+        variable_costs : numeric (iterable or scalar)
+            The costs associated with one unit of the flow. If this is set the
+            costs will be added to the objective expression of the optimization
+            problem.
+        positive_gradient : :obj:`dict`, default: `{'ub': None, 'costs': 0}`
+            A dictionary containing the following two keys:
+
+             * `'ub'`: numeric (iterable, scalar or None), the normed *upper
+               bound* on the positive difference (`flow[t-1] < flow[t]`) of
+               two consecutive flow values.
+             * `'costs``: numeric (scalar or None), the gradient cost per
+               unit.
+
+        negative_gradient : :obj:`dict`, default: `{'ub': None, 'costs': 0}`
+            A dictionary containing the following two keys:
+
+              * `'ub'`: numeric (iterable, scalar or None), the normed *upper
+                bound* on the negative difference (`flow[t-1] > flow[t]`) of
+                two consecutive flow values.
+              * `'costs``: numeric (scalar or None), the gradient cost per
+                unit.
+        """
+
+        def __init__(self, **kwargs):
+            dictionaries = ["positive_gradient", "negative_gradient"]
+            defaults = {
+                "positive_gradient": {"ub": None, "costs": 0},
+                "negative_gradient": {"ub": None, "costs": 0},
+            }
+
+            self.variable_costs = sequence(kwargs.get("variable_costs", 0))
+
+            for attribute in dictionaries:
+                value = kwargs.get(attribute, defaults.get(attribute))
+
+                setattr(
+                    self,
+                    attribute,
+                    {"ub": sequence(value["ub"]), "costs": value["costs"]},
+                )
diff --git a/tests/test_scripts/test_solph/test_multiobjective/test_model.py b/tests/test_scripts/test_solph/test_multiobjective/test_model.py
new file mode 100644
index 000000000..f2a3f0682
--- /dev/null
+++ b/tests/test_scripts/test_solph/test_multiobjective/test_model.py
@@ -0,0 +1,121 @@
+# -*- coding: utf-8 -*-
+"""
+This script contains tests for the MultiObjectiveModel
+"""
+
+import pandas as pd
+from pyomo.core.expr import current
+
+import oemof.solph as solph
+from oemof.solph.options import MultiObjective as mo
+
+
+def test_objective_functions_keys():
+
+    t = pd.date_range(start="01/01/2020", end="01/02/2020", freq="12H")
+
+    es = solph.EnergySystem(timeindex=t)
+
+    bus = solph.Bus(label='bus')
+
+    flow1 = solph.Flow(
+        label='flow1',
+        multiobjective=mo(
+            eco=mo.Objective(
+                variable_costs=1),
+            cost=mo.Objective(
+                variable_costs=2)))
+    flow2 = solph.Flow(
+        label='flow2',
+        multiobjective=mo(
+            eco=mo.Objective(
+                variable_costs=1),
+            benefit=mo.Objective(
+                variable_costs=3)))
+
+    trans = solph.Transformer(label='trans',
+                              inputs={bus: flow1},
+                              outputs={bus: flow2})
+
+    es.add(trans, bus)
+
+    model = solph.MultiObjectiveModel(es)
+
+    objective_functions = model.objective_functions
+
+    dict_compare = {'cost': 1, 'eco': 2, 'benefit': 3, '_standard': 4}
+
+    assert objective_functions.keys() == dict_compare.keys()
+
+
+def test_objective_functions_values():
+
+    t = pd.date_range(start="01/01/2020", end="01/02/2020", freq="12H")
+
+    es = solph.EnergySystem(timeindex=t)
+
+    bus = solph.Bus(label='bus')
+
+    costs = {'eco': 1, 'cost': 2, 'benefit': 3, '_standard': 0}
+
+    flow1 = solph.Flow(
+        label='flow1',
+        multiobjective=mo(
+            eco=mo.Objective(
+                variable_costs=costs['eco']),
+            cost=mo.Objective(
+                variable_costs=costs['cost'])))
+    flow2 = solph.Flow(
+        label='flow2',
+        multiobjective=mo(
+            eco=mo.Objective(
+                variable_costs=costs['eco']),
+            benefit=mo.Objective(
+                variable_costs=costs['benefit'])))
+
+    trans = solph.Transformer(
+        label='trans',
+        inputs={bus: flow1},
+        outputs={bus: flow2})
+
+    es.add(trans, bus)
+
+    model = solph.MultiObjectiveModel(es)
+
+    objective_functions = model.objective_functions
+
+    expr_compare = {
+        'eco': ('12.0*flow[trans,bus,0] + 12.0*flow[trans,bus,1] + '
+                + '12.0*flow[trans,bus,2] + 12.0*flow[bus,trans,0] + '
+                + '12.0*flow[bus,trans,1] + 12.0*flow[bus,trans,2]'),
+        'cost': ('24.0*flow[bus,trans,0] + 24.0*flow[bus,trans,1] + '
+                 + '24.0*flow[bus,trans,2]'),
+        'benefit': ('36.0*flow[trans,bus,0] + 36.0*flow[trans,bus,1] + '
+                    + '36.0*flow[trans,bus,2]'),
+        '_standard': '0'}
+
+    for x, c in costs.items():
+        mod_expr = current.expression_to_string(objective_functions[x])
+
+        assert mod_expr == expr_compare[x]
+
+
+def test_objective_functions_key_standard():
+
+    t = pd.date_range(start="01/01/2020", end="01/02/2020", freq="12H")
+
+    es = solph.EnergySystem(timeindex=t)
+
+    bus = solph.Bus(label='bus')
+
+    snk = solph.Sink(label='snk', inputs={bus: solph.Flow()})
+
+    es.add(snk, bus)
+
+    model = solph.MultiObjectiveModel(es)
+
+    objective_functions = model.objective_functions
+
+    dict_compare = {'_standard': 0}
+
+    assert objective_functions.keys() == dict_compare.keys()
diff --git a/tests/test_scripts/test_solph/test_multiobjective/test_model_solve.py b/tests/test_scripts/test_solph/test_multiobjective/test_model_solve.py
new file mode 100644
index 000000000..ff35fbd53
--- /dev/null
+++ b/tests/test_scripts/test_solph/test_multiobjective/test_model_solve.py
@@ -0,0 +1,107 @@
+# -*- coding: utf-8 -*-
+"""
+This script cains tests for the MultiObjectiveModel solve-function.
+"""
+
+import pandas as pd
+from nose.tools import raises
+from pyomo.core.expr import current
+
+import oemof.solph as solph
+from oemof.solph.options import MultiObjective as mo
+
+
+def test_opt_type():
+    es = solph.EnergySystem(timeindex=pd.date_range(
+        start='1/1/2018', end='3/1/2018', freq='D'))
+    mod = solph.MultiObjectiveModel(es)
+    with raises(Exception, match="Invalid optimization type"):
+        mod.solve(solver="cbc", optimization_type="flux_capacitator")
+
+
+def test_singular_obj_type():
+    es = solph.EnergySystem(timeindex=pd.date_range(
+        start='1/1/2018', end='3/1/2018', freq='D'))
+    mod = solph.MultiObjectiveModel(es)
+    with raises(TypeError, match='Objective is not of type "string"'):
+        mod.solve(solver="cbc",
+                  optimization_type='singular',
+                  objective=['eco'])
+
+
+def test_weighted_weights_dict():
+    es = solph.EnergySystem(timeindex=pd.date_range(
+        start='1/1/2018', end='3/1/2018', freq='D'))
+    bus = solph.Bus(label='bus')
+    src = solph.Source(label='src', outputs={bus: solph.Flow(
+        multiobjective=mo(
+            eco=mo.Objective(
+                variable_costs=0.5)))})
+    es.add(bus)
+    es.add(src)
+    mod = solph.MultiObjectiveModel(es)
+    with raises(TypeError, match="Objective weights must be of type 'dict'"):
+        mod.solve(solver="cbc",
+                  optimization_type='weighted',
+                  objective_weights=[0.3])
+
+
+def test_weighted_weights_length():
+    es = solph.EnergySystem(timeindex=pd.date_range(
+        start='1/1/2018', end='3/1/2018', freq='D'))
+    bus = solph.Bus(label='bus')
+    src = solph.Source(label='src', outputs={bus: solph.Flow(
+        multiobjective=mo(
+            eco=mo.Objective(
+                variable_costs=0.5)))})
+    es.add(bus)
+    es.add(src)
+    mod = solph.MultiObjectiveModel(es)
+    with raises(ValueError, match='Objective weights must not be empty'):
+        mod.solve(solver="cbc",
+                  optimization_type='weighted',
+                  objective_weights={})
+
+
+def test_objective_weighting():
+    es = solph.EnergySystem(timeindex=pd.date_range(
+        start='1/1/2018', end='3/1/2018', freq='D'))
+    bus = solph.Bus(label='bus')
+    src = solph.Source(label='src', outputs={bus: solph.Flow(
+        multiobjective=mo(
+            eco=mo.Objective(
+                variable_costs=0.5),
+            mon=mo.Objective(
+                variable_costs=1.3)))})
+    snk = solph.Sink(label='snk', inputs={bus: solph.Flow()})
+    es.add(bus)
+    es.add(src)
+    es.add(snk)
+    mod = solph.MultiObjectiveModel(es)
+    mod.solve(solver='cbc',
+              optimization_type='weighted',
+              objective_weights={'eco': 0.3, 'mon': 0.7})
+    s = current.expression_to_string(mod.objective)
+
+    # initialise expression string for first part of wighted sum
+    i = 0
+    expr_str = '(0.3*({0:}*flow[src,bus,{1:d}]'.format(0.5 * 24, i)
+    i += 1
+    # add summed values
+    while i < (len(es.timeindex) - 1):
+        expr_str += ' + {0:}*flow[src,bus,{1:d}]'.format(0.5 * 24, i)
+        i += 1
+    # add end of first partial sum
+    expr_str += ' + {0:}*flow[src,bus,{1:d}])'.format(0.5 * 24, i)
+    # initialise expression string for second part of wighted sum
+    i = 0
+    expr_str += ' + 0.7*({0:}*flow[src,bus,{1:d}]'.format(1.3 * 24, i)
+    i += 1
+    # add summed values
+    while i < (len(es.timeindex) - 1):
+        expr_str += ' + {0:}*flow[src,bus,{1:d}]'.format(1.3 * 24, i)
+        i += 1
+    # add end of first partial sum
+    expr_str += ' + {0:}*flow[src,bus,{1:d}]))'.format(1.3 * 24, i)
+
+    assert s == expr_str
diff --git a/tests/test_scripts/test_solph/test_multiobjective/test_multiobjective.py b/tests/test_scripts/test_solph/test_multiobjective/test_multiobjective.py
new file mode 100644
index 000000000..afac70cac
--- /dev/null
+++ b/tests/test_scripts/test_solph/test_multiobjective/test_multiobjective.py
@@ -0,0 +1,140 @@
+# -*- coding: utf-8 -*-
+"""
+This script contains an integration test for multi objective modelling.
+"""
+
+import pandas as pd
+
+import oemof.solph as solph
+from oemof.solph import processing
+from oemof.solph import views
+from oemof.solph.options import MultiObjective as mo
+
+
+def test_multiobjective():
+    ###########################################################################
+    # Set Costs, Demand etc.
+    ###########################################################################
+
+    # set list with different weights
+    weight_list = [(1, 0), (0.75, 0.25), (0.5, 0.5), (0.25, 0.75), (0, 1)]
+
+    # set fixed costs
+    cost_ecological_fix = 2
+    cost_financial_fix = 2
+
+    # read import data
+    import_data = pd.DataFrame(
+        index=pd.date_range(start='2020-01-01', periods=8, freq='3H'),
+        data={'cost_ecological': [0, 1, 1, 1, 1, 2, 3, 4],
+              'cost_financial': [0, 1, 2, 3, 4, 1, 1, 1],
+              'p_el_dem': [1, 1, 1, 1, 1, 1, 1, 1]})
+
+    # get demand
+    p_el_dem = import_data.loc[:, 'p_el_dem']
+
+    # get costs
+    cost_ecological_var = import_data.loc[:, 'cost_ecological']
+    cost_financial_var = import_data.loc[:, 'cost_financial']
+
+    ###########################################################################
+    # Create Enegy System
+    ###########################################################################
+    # get timeindex
+    timesteps = pd.date_range(
+        start=import_data.index[0],
+        end=import_data.index[-1],
+        freq=import_data.index.inferred_freq)
+
+    # initialise energy system
+    energy_system = solph.EnergySystem(timeindex=timesteps)
+
+    # create and add electrical bus
+    el_bus = solph.Bus(
+        label='el_bus')
+    energy_system.add(el_bus)
+
+    # create and add electrical demand
+    el_dem = solph.Sink(
+        label='el_dem',
+        inputs={el_bus: solph.Flow(fix=p_el_dem, nominal_value=1000)})
+    energy_system.add(el_dem)
+
+    # create and add electrical source with variable price
+    el_source_var = solph.Source(
+        label='el_source_var',
+        outputs={el_bus: solph.Flow(multiobjective=mo(
+            ecological=mo.Objective(
+                variable_costs=cost_ecological_var),
+            financial=mo.Objective(
+                variable_costs=cost_financial_var)))})
+    energy_system.add(el_source_var)
+
+    # create and add electrical source with fixed price
+    el_source_fix = solph.Source(
+        label='el_source_fix',
+        outputs={el_bus: solph.Flow(multiobjective=mo(
+            ecological=mo.Objective(
+                variable_costs=cost_ecological_fix),
+            financial=mo.Objective(
+                variable_costs=cost_financial_fix)))})
+    energy_system.add(el_source_fix)
+
+    #######################################################################
+    # Create LP Model
+    #######################################################################
+
+    # create a pyomo optimization problem
+    optimisation_model = solph.MultiObjectiveModel(energy_system)
+
+    ###########################################################################
+    # Iterate over different weights and Solve LP Model
+    ###########################################################################
+
+    for weight_ecological, weight_financial in weight_list:
+
+        # solve problem using cplex
+        # no console output
+        optimisation_model.solve(
+            solver='cbc',
+            optimization_type='weighted',
+            objective_weights={'ecological': weight_ecological,
+                               'financial': weight_financial},
+            solve_kwargs={'tee': False},
+            cmdline_options={'tmlim': 30})
+
+        #######################################################################
+        # Extract and Postprocess results
+        #######################################################################
+
+        # get result dict
+        results = processing.results(optimisation_model)
+
+        # get results for important nodes
+        el_bus_data = views.node(results, 'el_bus')['sequences']
+
+        # get costs for variable cost source
+        p_el_var = el_bus_data[(('el_source_var', 'el_bus'), 'flow')]
+
+        # get costs for fixed cost source
+        p_el_fix = el_bus_data[(('el_source_fix', 'el_bus'), 'flow')]
+
+        # get variable costs
+        cost_var = (weight_ecological * cost_ecological_var
+                    + weight_financial * cost_financial_var)
+
+        # get fix costs
+        cost_fix = (weight_ecological * cost_ecological_fix
+                    + weight_financial * cost_financial_fix)
+
+        ###################################################################
+        # Validate results
+        ###################################################################
+
+        # iterate over timesteps
+        for t in range(len(timesteps)):
+            # assert correct assignment only for unambiguous timesteps
+            if cost_var.iat[t] > cost_fix:
+                assert p_el_var.iat[t] == 0
+            elif cost_var.iat[t] < cost_fix:
+                assert p_el_fix.iat[t] == 0