diff --git a/examples/knowledgebases/__init__.py b/examples/knowledgebases/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/examples/knowledgebases/ssbj/Cnstrnt_sigmas.py b/examples/knowledgebases/ssbj/Cnstrnt_sigmas.py
new file mode 100644
index 0000000000000000000000000000000000000000..63e8e24f62d6b876acbca8cb99bd2b6661c228ad
--- /dev/null
+++ b/examples/knowledgebases/ssbj/Cnstrnt_sigmas.py
@@ -0,0 +1,116 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+SSBJ test case - http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19980234657.pdf
+Original Python implementation for OpenMDAO integration developed by
+Sylvain Dubreuil and Remi Lafage of ONERA, the French Aerospace Lab.
+Orignal files taken from: https://github.com/OneraHub/SSBJ-OpenMDAO
+The files were adjusted for optimal use in KADMOS by Imco van Gent.
+"""
+from __future__ import absolute_import, division, print_function
+
+import numpy as np
+from lxml import etree
+
+from examples.knowledgebases.ssbj import root_tag, x_tc, x_AR, x_Sref, x_lambda, x_section, x_L, \
+    x_sigma1, x_sigma2, x_sigma3, x_sigma4, x_sigma5
+from examples.knowledgebases.ssbj.common import PolynomialFunction, add_discipline_to_cmdows
+from openlego.api import AbstractDiscipline
+from openlego.utils.xml_utils import xml_safe_create_element
+
+
+class Cnstrnt_sigmas(AbstractDiscipline):
+
+    @property
+    def creator(self):
+        return u'S. Dubreuil and R. Lafage'
+
+    @property
+    def owner(self):
+        return u'J. Sobieszczanski-Sobieski'
+
+    @property
+    def operator(self):
+        return u'I. van Gent'
+
+    @property
+    def description(self):
+        return u'Structural stress constraint of the SSBJ test case.'
+
+    @property
+    def description(self):
+        return u'First discipline of the Sellar problem'
+
+    def generate_input_xml(self):
+        root = etree.Element(root_tag)
+        doc = etree.ElementTree(root)
+
+        xml_safe_create_element(doc, x_tc, 0.05)
+        xml_safe_create_element(doc, x_AR, 5.5)
+        xml_safe_create_element(doc, x_Sref, 1000.0)
+        xml_safe_create_element(doc, x_lambda, 0.25)
+        xml_safe_create_element(doc, x_section, 1.0)
+        xml_safe_create_element(doc, x_L, 49909.58578)
+
+        return etree.tostring(doc, encoding='utf-8', pretty_print=True, xml_declaration=True)
+
+    def generate_output_xml(self):
+        root = etree.Element(root_tag)
+        doc = etree.ElementTree(root)
+
+        xml_safe_create_element(doc, x_sigma1, 1.12255)
+        xml_safe_create_element(doc, x_sigma2, 1.08170213)
+        xml_safe_create_element(doc, x_sigma3, 1.0612766)
+        xml_safe_create_element(doc, x_sigma4, 1.04902128)
+        xml_safe_create_element(doc, x_sigma5, 1.04085106)
+
+        return etree.tostring(doc, encoding='utf-8', pretty_print=True, xml_declaration=True)
+
+    def deploy(self):
+        """Deploy this discipline's template in-/output XML files and its information in the CMDOWS file."""
+        with open(self.in_file, 'w') as f:
+            f.write(self.generate_input_xml())
+        with open(self.out_file, 'w') as f:
+            f.write(self.generate_output_xml())
+        add_discipline_to_cmdows(self)
+
+    @staticmethod
+    def execute(in_file, out_file):
+        doc = etree.parse(in_file)
+        z0 = float(doc.xpath(x_tc)[0].text)
+        z3 = float(doc.xpath(x_AR)[0].text)
+        z5 = float(doc.xpath(x_Sref)[0].text)
+        x0 = float(doc.xpath(x_lambda)[0].text)
+        x1 = float(doc.xpath(x_section)[0].text)
+        L = float(doc.xpath(x_L)[0].text)
+
+        pf = PolynomialFunction()
+
+        b = np.sqrt(abs(z5 * z3)) / 2.0
+        R = (1.0 + 2.0 * x0) / (3.0 * (1.0 + x0))
+
+        Sigma0 = pf.eval([z0, L, x1, b, R], [4, 1, 4, 1, 1], [0.1] * 5, "sigma[1]")
+        Sigma1 = pf.eval([z0, L, x1, b, R], [4, 1, 4, 1, 1], [0.15] * 5, "sigma[2]")
+        Sigma2 = pf.eval([z0, L, x1, b, R], [4, 1, 4, 1, 1], [0.2] * 5, "sigma[3]")
+        Sigma3 = pf.eval([z0, L, x1, b, R], [4, 1, 4, 1, 1], [0.25] * 5, "sigma[4]")
+        Sigma4 = pf.eval([z0, L, x1, b, R], [4, 1, 4, 1, 1], [0.30] * 5, "sigma[5]")
+
+        root = etree.Element(root_tag)
+        doc = etree.ElementTree(root)
+        xml_safe_create_element(doc, x_sigma1, Sigma0)
+        xml_safe_create_element(doc, x_sigma2, Sigma1)
+        xml_safe_create_element(doc, x_sigma3, Sigma2)
+        xml_safe_create_element(doc, x_sigma4, Sigma3)
+        xml_safe_create_element(doc, x_sigma5, Sigma4)
+        doc.write(out_file, encoding='utf-8', pretty_print=True, xml_declaration=True)
+
+
+if __name__ == "__main__":
+
+    cnstrnt_sigmas_analysis = Cnstrnt_sigmas()
+
+    in_file = '__test__cnstrnt_sigmas_input.xml'
+    out_file = '__test__cnstrnt_sigmas_output.xml'
+    with open(in_file, 'w') as f:
+        f.write(cnstrnt_sigmas_analysis.generate_input_xml())
+    cnstrnt_sigmas_analysis.execute(in_file, out_file)
diff --git a/examples/knowledgebases/ssbj/Cnstrnt_theta.py b/examples/knowledgebases/ssbj/Cnstrnt_theta.py
new file mode 100644
index 0000000000000000000000000000000000000000..e46e5dfd6ad402167f44ff3c9f9a7c19ac489480
--- /dev/null
+++ b/examples/knowledgebases/ssbj/Cnstrnt_theta.py
@@ -0,0 +1,96 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+SSBJ test case - http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19980234657.pdf
+Original Python implementation for OpenMDAO integration developed by
+Sylvain Dubreuil and Remi Lafage of ONERA, the French Aerospace Lab.
+Orignal files taken from: https://github.com/OneraHub/SSBJ-OpenMDAO
+The files were adjusted for optimal use in KADMOS by Imco van Gent.
+"""
+from __future__ import absolute_import, division, print_function
+
+import numpy as np
+from lxml import etree
+
+from examples.knowledgebases.ssbj import root_tag, x_AR, x_Sref, x_lambda, x_section, x_L, x_Theta
+from examples.knowledgebases.ssbj.common import PolynomialFunction, add_discipline_to_cmdows
+from openlego.api import AbstractDiscipline
+from openlego.utils.xml_utils import xml_safe_create_element
+
+
+class Cnstrnt_theta(AbstractDiscipline):
+
+    @property
+    def creator(self):
+        return u'S. Dubreuil and R. Lafage'
+
+    @property
+    def owner(self):
+        return u'J. Sobieszczanski-Sobieski'
+
+    @property
+    def operator(self):
+        return u'I. van Gent'
+
+    @property
+    def description(self):
+        return u'Structural constraint of the SSBJ test case.'
+
+    def generate_input_xml(self):
+        root = etree.Element(root_tag)
+        doc = etree.ElementTree(root)
+
+        xml_safe_create_element(doc, x_AR, 5.5)
+        xml_safe_create_element(doc, x_Sref, 1000.0)
+        xml_safe_create_element(doc, x_lambda, 0.25)
+        xml_safe_create_element(doc, x_section, 1.0)
+        xml_safe_create_element(doc, x_L, 49909.58578)
+
+        return etree.tostring(doc, encoding='utf-8', pretty_print=True, xml_declaration=True)
+
+    def generate_output_xml(self):
+        root = etree.Element(root_tag)
+        doc = etree.ElementTree(root)
+
+        xml_safe_create_element(doc, x_Theta, 0.950978)
+
+        return etree.tostring(doc, encoding='utf-8', pretty_print=True, xml_declaration=True)
+
+    def deploy(self):
+        """Deploy this discipline's template in-/output XML files and its information in the CMDOWS file."""
+        with open(self.in_file, 'w') as f:
+            f.write(self.generate_input_xml())
+        with open(self.out_file, 'w') as f:
+            f.write(self.generate_output_xml())
+        add_discipline_to_cmdows(self)
+
+    @staticmethod
+    def execute(in_file, out_file):
+        doc = etree.parse(in_file)
+        z3 = float(doc.xpath(x_AR)[0].text)
+        z5 = float(doc.xpath(x_Sref)[0].text)
+        x0 = float(doc.xpath(x_lambda)[0].text)
+        x1 = float(doc.xpath(x_section)[0].text)
+        L = float(doc.xpath(x_L)[0].text)
+
+        pf = PolynomialFunction()
+
+        b = np.sqrt(abs(z5 * z3)) / 2.0
+        R = (1.0 + 2.0 * x0) / (3.0 * (1.0 + x0))
+        Theta = pf.eval([abs(x1), b, R, L], [2, 4, 4, 3], [0.25] * 4, "twist")
+
+        root = etree.Element(root_tag)
+        doc = etree.ElementTree(root)
+        xml_safe_create_element(doc, x_Theta, Theta)
+        doc.write(out_file, encoding='utf-8', pretty_print=True, xml_declaration=True)
+
+
+if __name__ == "__main__":
+
+    cnstrnt_theta = Cnstrnt_theta()
+
+    in_file = '__test__cnstrnt_theta_input.xml'
+    out_file = '__test__cnstrnt_theta_output.xml'
+    with open(in_file, 'w') as f:
+        f.write(cnstrnt_theta.generate_input_xml())
+    cnstrnt_theta.execute(in_file, out_file)
diff --git a/examples/knowledgebases/ssbj/Structures.py b/examples/knowledgebases/ssbj/Structures.py
new file mode 100644
index 0000000000000000000000000000000000000000..3f5ef91e0472cf6b19f4c959c953efdcd13d31e4
--- /dev/null
+++ b/examples/knowledgebases/ssbj/Structures.py
@@ -0,0 +1,120 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+SSBJ test case - http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19980234657.pdf
+Original Python implementation for OpenMDAO integration developed by
+Sylvain Dubreuil and Remi Lafage of ONERA, the French Aerospace Lab.
+Orignal files taken from: https://github.com/OneraHub/SSBJ-OpenMDAO
+The files were adjusted for optimal use in KADMOS by Imco van Gent.
+"""
+from __future__ import absolute_import, division, print_function
+
+import numpy as np
+from lxml import etree
+
+from examples.knowledgebases.ssbj import root_tag, x_tc, x_AR, x_Lambda, x_Sref, x_lambda, x_section, x_WO, x_WE, x_WFO, \
+    x_L, x_Nz, x_WT, x_WF
+from examples.knowledgebases.ssbj.common import PolynomialFunction, add_discipline_to_cmdows
+from openlego.api import AbstractDiscipline
+from openlego.utils.xml_utils import xml_safe_create_element
+
+
+class Structures(AbstractDiscipline):
+
+    @property
+    def creator(self):
+        return u'S. Dubreuil and R. Lafage'
+
+    @property
+    def owner(self):
+        return u'J. Sobieszczanski-Sobieski'
+
+    @property
+    def operator(self):
+        return u'I. van Gent'
+
+    @property
+    def description(self):
+        return u'Structural analysis discipline of the SSBJ test case.'
+
+    def generate_input_xml(self):
+        root = etree.Element(root_tag)
+        doc = etree.ElementTree(root)
+
+        xml_safe_create_element(doc, x_tc, 0.05)
+        xml_safe_create_element(doc, x_AR, 5.5)
+        xml_safe_create_element(doc, x_Lambda, 55.0)
+        xml_safe_create_element(doc, x_Sref, 1000.0)
+        xml_safe_create_element(doc, x_lambda, 0.25)
+        xml_safe_create_element(doc, x_section, 1.0)
+        xml_safe_create_element(doc, x_WO, 25000.)
+        xml_safe_create_element(doc, x_WE, 5748.915355)
+        xml_safe_create_element(doc, x_WFO, 2000.)
+        xml_safe_create_element(doc, x_L, 49909.58578)
+        xml_safe_create_element(doc, x_Nz, 6.0)
+
+        return etree.tostring(doc, encoding='utf-8', pretty_print=True, xml_declaration=True)
+
+    def generate_output_xml(self):
+        root = etree.Element(root_tag)
+        doc = etree.ElementTree(root)
+
+        xml_safe_create_element(doc, x_WT, 49909.58578)
+        xml_safe_create_element(doc, x_WF, 7306.20261)
+
+        return etree.tostring(doc, encoding='utf-8', pretty_print=True, xml_declaration=True)
+
+    def deploy(self):
+        """Deploy this discipline's template in-/output XML files and its information in the CMDOWS file."""
+        with open(self.in_file, 'w') as f:
+            f.write(self.generate_input_xml())
+        with open(self.out_file, 'w') as f:
+            f.write(self.generate_output_xml())
+        add_discipline_to_cmdows(self)
+
+    @staticmethod
+    def execute(in_file, out_file):
+        doc = etree.parse(in_file)
+        z0 = float(doc.xpath(x_tc)[0].text)
+        z3 = float(doc.xpath(x_AR)[0].text)
+        z4 = float(doc.xpath(x_Lambda)[0].text)
+        z5 = float(doc.xpath(x_Sref)[0].text)
+        x0 = float(doc.xpath(x_lambda)[0].text)
+        x1 = float(doc.xpath(x_section)[0].text)
+        L = float(doc.xpath(x_L)[0].text)
+        WE = float(doc.xpath(x_WE)[0].text)
+        NZ = float(doc.xpath(x_Nz)[0].text)
+        WO = float(doc.xpath(x_WO)[0].text)
+        WFO = float(doc.xpath(x_WFO)[0].text)
+
+        pf = PolynomialFunction()
+
+        t = z0 * z5 / (np.sqrt(abs(z5 * z3)))
+
+        Fo1 = pf.eval([x1], [1], [.008], "Fo1")
+
+        WT_hat = L
+        WW = Fo1 * (0.0051 * abs(WT_hat * NZ) ** 0.557 * \
+                    abs(z5) ** 0.649 * abs(z3) ** 0.5 * abs(z0) ** (-0.4) \
+                    * abs(1.0 + x0) ** 0.1 * (0.1875 * abs(z5)) ** 0.1 \
+                    / abs(np.cos(z4 * np.pi / 180.)))
+        WFW = 5.0 / 18.0 * abs(z5) * 2.0 / 3.0 * t * 42.5
+        WF = WFW + WFO
+        WT = WO + WW + WF + WE
+
+        root = etree.Element(root_tag)
+        doc = etree.ElementTree(root)
+        xml_safe_create_element(doc, x_WF, WF)
+        xml_safe_create_element(doc, x_WT, WT)
+        doc.write(out_file, encoding='utf-8', pretty_print=True, xml_declaration=True)
+
+
+if __name__ == "__main__":
+
+    structural_analysis = Structures()
+
+    in_file = '__test__structures_core_input.xml'
+    out_file = '__test__structures_core_output.xml'
+    with open(in_file, 'w') as f:
+        f.write(structural_analysis.generate_input_xml())
+    structural_analysis.execute(in_file, out_file)
diff --git a/examples/knowledgebases/ssbj/__init__.py b/examples/knowledgebases/ssbj/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..58679895cf833f01652b0b93a0732c5f4ab0d391
--- /dev/null
+++ b/examples/knowledgebases/ssbj/__init__.py
@@ -0,0 +1,121 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+Copyright 2017 D. de Vries
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+
+This file contains code to clean and deploy the knowledge base of the test Sellar case.
+"""
+from __future__ import absolute_import, division, print_function
+
+import os
+import sys
+from shutil import copyfile
+
+from openlego.utils.xml_utils import xml_merge
+
+dir_path = os.path.dirname(os.path.realpath(__file__))
+base_file_path = os.path.join(dir_path, 'sellar-base.xml')
+
+root_tag = 'data_schema'
+cat1 = 'aircraft/geometry'
+cat2 = 'aircraft/weight'
+cat3 = 'aircraft/other'
+cat4 = 'aircraft/reference'
+x_root = '/' + root_tag
+
+x_tc = '/'.join([x_root, cat1, 'tc'])
+x_AR = '/'.join([x_root, cat1, 'AR'])
+x_Lambda = '/'.join([x_root, cat1, 'Lambda'])
+x_Sref = '/'.join([x_root, cat1, 'Sref'])
+x_Theta = '/'.join([x_root, cat1, 'Theta'])
+x_lambda = '/'.join([x_root, cat1, 'lambda'])
+x_section = '/'.join([x_root, cat1, 'section'])
+
+x_WT = '/'.join([x_root, cat2, 'WT'])
+x_WBE = '/'.join([x_root, cat2, 'WBE'])
+x_WE = '/'.join([x_root, cat2, 'WE'])
+x_WF = '/'.join([x_root, cat2, 'WF'])
+x_WFO = '/'.join([x_root, cat2, 'WFO'])
+x_WO = '/'.join([x_root, cat2, 'WO'])
+
+x_D = '/'.join([x_root, cat3, 'D'])
+x_L = '/'.join([x_root, cat3, 'L'])
+x_Cf = '/'.join([x_root, cat3, 'Cf'])
+x_CDmin = '/'.join([x_root, cat3, 'CDmin'])
+x_T = '/'.join([x_root, cat3, 'T'])
+x_DT = '/'.join([x_root, cat3, 'DT'])
+x_fin = '/'.join([x_root, cat3, 'fin'])
+x_SFC = '/'.join([x_root, cat3, 'SFC'])
+x_dpdx = '/'.join([x_root, cat3, 'dpdx'])
+x_R = '/'.join([x_root, cat3, 'R'])
+x_Nz = '/'.join([x_root, cat3, 'Nz'])
+x_sigma1 = '/'.join([x_root, cat3, 'sigma1'])
+x_sigma2 = '/'.join([x_root, cat3, 'sigma2'])
+x_sigma3 = '/'.join([x_root, cat3, 'sigma3'])
+x_sigma4 = '/'.join([x_root, cat3, 'sigma4'])
+x_sigma5 = '/'.join([x_root, cat3, 'sigma5'])
+
+x_h = '/'.join([x_root, cat4, 'h'])
+x_M = '/'.join([x_root, cat4, 'M'])
+x_ESF = '/'.join([x_root, cat4, 'ESF'])
+x_Temp = '/'.join([x_root, cat4, 'Temp'])
+
+from examples.knowledgebases.ssbj.Cnstrnt_sigmas import Cnstrnt_sigmas
+from examples.knowledgebases.ssbj.Cnstrnt_theta import Cnstrnt_theta
+from examples.knowledgebases.ssbj.Structures import Structures
+
+
+def list_disciplines():
+    return [Structures(), Cnstrnt_theta(), Cnstrnt_sigmas()]
+
+def try_to_remove(file):
+    try:
+        os.remove(file)
+    except:
+        pass
+
+def clean():
+    for discipline in list_disciplines():
+        try_to_remove(discipline.in_file)
+        try_to_remove(discipline.out_file)
+        try_to_remove(discipline.json_file)
+    try_to_remove(base_file_path)
+
+    dir = os.path.dirname(os.path.realpath(__file__))
+    for file in os.listdir(dir):
+        if '__test__' in file:
+            os.remove(file)
+        if '__cmdows__' in file:
+            os.remove(file)
+
+
+def deploy():
+    _copy = True
+    for discipline in list_disciplines():
+        discipline.deploy()
+        if _copy:
+            _copy = False
+            copyfile(discipline.in_file, base_file_path)
+        else:
+            xml_merge(base_file_path, discipline.in_file)
+
+        xml_merge(base_file_path, discipline.out_file)
+
+
+if __name__ == '__main__':
+    if len(sys.argv) == 1:
+        deploy()
+    elif len(sys.argv) == 2 and sys.argv[1] == 'clean':
+        clean()
diff --git a/examples/knowledgebases/ssbj/common.py b/examples/knowledgebases/ssbj/common.py
new file mode 100644
index 0000000000000000000000000000000000000000..5de52730be4a73bcdaf6859716e7d8dd4cbd10c0
--- /dev/null
+++ b/examples/knowledgebases/ssbj/common.py
@@ -0,0 +1,124 @@
+"""
+SSBJ test case - http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19980234657.pdf
+Python implementation and OpenMDAO integration developed by
+Sylvain Dubreuil and Remi Lafage of ONERA, the French Aerospace Lab.
+"""
+import os
+import numpy as np
+# pylint: disable=C0103
+from kadmos.cmdows import CMDOWS
+
+
+# TODO: This function should be part of OpenLEGO
+def add_discipline_to_cmdows(Discipline):
+    # Database files
+    files = os.listdir(os.path.dirname(os.path.realpath(__file__)))
+    cmdows_files = [file for file in files if '__cmdows__' in file]
+    assert len(cmdows_files) <= 1, 'Multiple CMDOWS files were found {}.'.format(cmdows_files)
+    if cmdows_files:
+        cmdows = CMDOWS(cmdows_files[0])
+    else:
+        cmdows_files = ['__cmdows__SSBJ.xml']
+        cmdows = CMDOWS()
+        cmdows.add_header()
+    cmdows.add_dc(Discipline.name, Discipline.name, 'main', 1, Discipline.version, Discipline.name)
+    cmdows.save(cmdows_files[0], pretty_print=True)
+
+
+class PolynomialFunction(object):
+
+    R = [[0.2736, 0.3970, 0.8152, 0.9230, 0.1108],
+         [0.4252, 0.4415, 0.6357, 0.7435, 0.1138],
+         [0.0329, 0.8856, 0.8390, 0.3657, 0.0019],
+         [0.0878, 0.7248, 0.1978, 0.0200, 0.0169],
+         [0.8955, 0.4568, 0.8075, 0.9239, 0.2525]]
+
+    d = dict()
+
+    def eval(self, S_new, flag, S_bound, var, deriv=False):
+        if len(S_new) > 1:
+            res = np.array([])
+            S_new = np.append(res, S_new)
+
+        if var not in self.d:
+            self.d[var] = list(S_new)
+
+        S = self.d[var]
+        S_norm = []
+        S_shifted = []
+        Ai = []
+        Aij = [[0.0]*len(S_new) for i in range(len(S_new))]
+
+        assert len(S) == len(S_new)
+
+        for i in range(len(S)):
+            S_norm.append(S_new[i] / S[i])
+
+            if S_norm[i] > 1.25:
+                S_norm[i] = 1.25
+            elif S_norm[i] < 0.75:
+                S_norm[i] = 0.75
+
+            S_shifted.append(S_norm[i]-1)
+
+            a = 0.1
+            b = a
+
+            if flag[i] == 3:
+                a = -a
+                b = a
+            elif flag[i] == 2:
+                b = 2*a
+            elif flag[i] == 4:
+                a = -a
+                b = 2*a
+
+            So = 0.0
+            Sl = So - S_bound[i]
+            Su = So + S_bound[i]
+            Mtx_shifted = np.array([[1.0, Sl, Sl**2],
+                                    [1.0, So, So**2],
+                                    [1.0, Su, Su**2]])
+
+            if flag[i] == 5:
+                F_bound = np.array([[1+(0.5*a)**2], [1.0], [1+(0.5*b)**2]])
+            else:
+                F_bound = np.array([[1-(0.5*a)], [1.0], [1+(0.5*b)]])
+
+            A = np.linalg.solve(Mtx_shifted, F_bound)
+
+            Ao = A[0]
+            B = A[1]
+
+            if var == "Fo1":
+                Ai.append(B)
+            else:
+                Ai.append(A[1])
+
+            Aij[i][i] = A[2]
+
+        for i in range(len(S)):
+            for j in range(i+1, len(S)):
+                Aij[i][j] = Aij[i][i] * self.R[i][j]
+                Aij[j][i] = Aij[i][j]
+
+        Ai = np.matrix(np.array(Ai))
+        Aij = np.matrix(np.array(Aij))
+        S_shifted = np.matrix(np.array(S_shifted))
+
+        if deriv:
+            return S_shifted, Ai, Aij
+        else:
+            return float((Ao + Ai.T * S_shifted.T + 0.5 * S_shifted * Aij * S_shifted.T)[0])
+
+
+if __name__ == '__main__':
+
+    p = PolynomialFunction()
+
+    init = [1.0, 37.080992435478315, 0.4, 1000.0]
+    b = [1.0, 37.080992435478315, 0.4, 26315.848165047268]
+    a = [1.0, 37.080992435478315, 0.4, -12243.514743699088]
+
+    print("it 1", p.eval([1.0], [1], [0.008], "Fo1"))
+    print("it 2", p.eval([0.766], [1], [0.008], "Fo1"))
diff --git a/examples/knowledgebases/ssbj/data-schema.xsd b/examples/knowledgebases/ssbj/data-schema.xsd
index fe4896065b6946c6c754a1af03b92cf9b55ec146..0c72e93dd25c247bb5ecc407f35674e2dfdc19c8 100644
--- a/examples/knowledgebases/ssbj/data-schema.xsd
+++ b/examples/knowledgebases/ssbj/data-schema.xsd
@@ -9,10 +9,14 @@
     <xs:element name="Lambda" type="xs:decimal"/>
     <xs:element name="Sref" type="xs:decimal"/>
     <xs:element name="Cf" type="xs:decimal"/>
+    <xs:element name="CDmin" type="xs:decimal"/>
     <xs:element name="T" type="xs:decimal"/>
     <xs:element name="WT" type="xs:decimal"/>
     <xs:element name="WF" type="xs:decimal"/>
     <xs:element name="WE" type="xs:decimal"/>
+    <xs:element name="WFO" type="xs:decimal"/>
+    <xs:element name="WBE" type="xs:decimal"/>
+    <xs:element name="WO" type="xs:decimal"/>
     <xs:element name="D" type="xs:decimal"/>
     <xs:element name="L" type="xs:decimal"/>
     <xs:element name="Theta" type="xs:decimal"/>
@@ -23,9 +27,14 @@
     <xs:element name="section" type="xs:decimal"/>
     <xs:element name="dpdx" type="xs:decimal"/>
     <xs:element name="R" type="xs:decimal"/>
+    <xs:element name="Nz" type="xs:decimal"/>
     <xs:element name="DT" type="xs:decimal"/>
     <xs:element name="Temp" type="xs:decimal"/>
-    <xs:element name="sigma" type="xs:string"/>
+    <xs:element name="sigma1" type="xs:decimal"/>
+    <xs:element name="sigma2" type="xs:decimal"/>
+    <xs:element name="sigma3" type="xs:decimal"/>
+    <xs:element name="sigma4" type="xs:decimal"/>
+    <xs:element name="sigma5" type="xs:decimal"/>
 
     <!-- definition of complex elements -->
     <xs:element name="geometry">
@@ -45,9 +54,12 @@
     <xs:element name="weight">
       <xs:complexType>
         <xs:all>
-            <xs:element ref="WE" minOccurs="0"/>
             <xs:element ref="WT" minOccurs="0"/>
+            <xs:element ref="WBE" minOccurs="0"/>
+            <xs:element ref="WE" minOccurs="0"/>
             <xs:element ref="WF" minOccurs="0"/>
+            <xs:element ref="WFO" minOccurs="0"/>
+            <xs:element ref="WO" minOccurs="0"/>
         </xs:all>
       </xs:complexType>
     </xs:element>
@@ -58,13 +70,19 @@
             <xs:element ref="D" minOccurs="0"/>
             <xs:element ref="L" minOccurs="0"/>
             <xs:element ref="Cf" minOccurs="0"/>
+            <xs:element ref="CDmin" minOccurs="0"/>
             <xs:element ref="T" minOccurs="0"/>
             <xs:element ref="DT" minOccurs="0"/>
             <xs:element ref="fin" minOccurs="0"/>
             <xs:element ref="SFC" minOccurs="0"/>
             <xs:element ref="dpdx" minOccurs="0"/>
             <xs:element ref="R" minOccurs="0"/>
-            <xs:element ref="sigma" minOccurs="0"/>
+            <xs:element ref="Nz" minOccurs="0"/>
+            <xs:element ref="sigma1" minOccurs="0"/>
+            <xs:element ref="sigma2" minOccurs="0"/>
+            <xs:element ref="sigma3" minOccurs="0"/>
+            <xs:element ref="sigma4" minOccurs="0"/>
+            <xs:element ref="sigma5" minOccurs="0"/>
         </xs:all>
       </xs:complexType>
     </xs:element>
diff --git a/examples/knowledgebases/ssbj/ssbj_toolrepo_cmdows.xml b/examples/knowledgebases/ssbj/ssbj_toolrepo_cmdolds.xml
similarity index 100%
rename from examples/knowledgebases/ssbj/ssbj_toolrepo_cmdows.xml
rename to examples/knowledgebases/ssbj/ssbj_toolrepo_cmdolds.xml
diff --git a/examples/knowledgebases/ssbj/structure-input.xml b/examples/knowledgebases/ssbj/structure-input.xml
deleted file mode 100644
index ce1d83d4d0f44d94b2a3fa2b449a44859623b209..0000000000000000000000000000000000000000
--- a/examples/knowledgebases/ssbj/structure-input.xml
+++ /dev/null
@@ -1,19 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<data_schema xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="data-schema.xsd">
-    <aircraft>
-        <geometry>
-            <tc>0.05</tc>
-            <AR>5.5</AR>
-            <Lambda>55.0</Lambda>
-            <Sref>1000.0</Sref>
-            <lambda>0.25</lambda>
-            <section>1.0</section>
-        </geometry>
-        <weight>
-            <WE>5748.915355</WE>
-        </weight>
-        <other>
-            <L>49909.58578</L>
-        </other>
-    </aircraft>
-</data_schema>
\ No newline at end of file
diff --git a/examples/knowledgebases/ssbj/structure-output.xml b/examples/knowledgebases/ssbj/structure-output.xml
deleted file mode 100644
index 692b725558a138e45f4dec60f981bb5693ed0430..0000000000000000000000000000000000000000
--- a/examples/knowledgebases/ssbj/structure-output.xml
+++ /dev/null
@@ -1,15 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<data_schema xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="data-schema.xsd">
-    <aircraft>
-        <geometry>
-            <Theta>0.950978</Theta>
-        </geometry>
-        <weight>
-            <WT>49909.58578</WT>
-            <WF>7306.20261</WF>
-        </weight>
-        <other>
-            <sigma>1.12255, 1.08170213, 1.0612766, 1.04902128, 1.04085106</sigma>
-        </other>
-    </aircraft>
-</data_schema>
\ No newline at end of file
diff --git a/examples/scripts/sellar_problem/(X)DSM/Mdao_CO.pdf b/examples/scripts/sellar_problem/(X)DSM/Mdao_CO.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..cbe115e25faa50bb9760d4c5d953032784b7d6b7
Binary files /dev/null and b/examples/scripts/sellar_problem/(X)DSM/Mdao_CO.pdf differ
diff --git a/examples/scripts/sellar_problem/(X)DSM/Test_XDSM.pdf b/examples/scripts/sellar_problem/(X)DSM/Test_XDSM.pdf
deleted file mode 100644
index 10f356bfb948d73001b4faa82e756750cef4cd4c..0000000000000000000000000000000000000000
Binary files a/examples/scripts/sellar_problem/(X)DSM/Test_XDSM.pdf and /dev/null differ
diff --git a/examples/scripts/sellar_problem/(X)DSM/Test_XDSM.tex b/examples/scripts/sellar_problem/(X)DSM/Test_XDSM.tex
deleted file mode 100644
index ed291652989a9a6e1ae35279abcf188a26072fd7..0000000000000000000000000000000000000000
--- a/examples/scripts/sellar_problem/(X)DSM/Test_XDSM.tex
+++ /dev/null
@@ -1,227 +0,0 @@
-\documentclass{article}
-\usepackage{geometry}
-\usepackage{amsfonts}
-\usepackage{amsmath}
-\usepackage{amssymb}
-\usepackage{tikz}
-
-\input{/Users/imcovangent/Documents/PhD/Software/KADMOS/kadmos/external/XDSM_writer/diagram_border}
-
-\begin{document}
-
-\input{/Users/imcovangent/Documents/PhD/Software/KADMOS/kadmos/external/XDSM_writer/diagram_styles}
-
-\begin{tikzpicture}
-
-  \matrix[MatrixSetup]
-  {
-    %Row 1
-    &
-    &
-    &
-    &
-    &
-    &
-    &
-    &
-    &
-    &
-    &
-    &
-    \\
-    &
-    \node [Coordinator] (436f6f7264696e61746f72) {COOR}; &
-    \node [DataIO] (41-436f6f7264696e61746f72) {$a$}; &
-    \node [DataIO] (53797374656d2d4f7074696d697a6572-436f6f7264696e61746f72) {$x1^{c,0}$\\[1pt] $y2^{c,0}$\\[1pt] $y1^{c,0}$\\[1pt] $z1^{0}$\\[1pt] $z2^{0}$}; &
-    &
-    \node [DataIO] (5375624f7074696d697a65722d30-436f6f7264696e61746f72) {$z2^{c,0}$\\[1pt] $z1^{c,0}$\\[1pt] $x1^{0}$}; &
-    &
-    &
-    &
-    \node [DataIO] (5375624f7074696d697a65722d31-436f6f7264696e61746f72) {$z1^{c,i2,0}$\\[1pt] $z2^{c,i2,0}$}; &
-    &
-    &
-    \\
-    &
-    &
-    \node [PreAnalysis] (41) {A}; &
-    &
-    &
-    &
-    \node [DataInter] (4431-41) {$c$}; &
-    &
-    &
-    &
-    \node [DataInter] (4432-41) {$c$}; &
-    &
-    \\
-    &
-    \node [DataIO] (436f6f7264696e61746f72-53797374656d2d4f7074696d697a6572) {$y1^{c,*}$\\[1pt] $x1^{c,*}$\\[1pt] $y2^{c,*}$\\[1pt] $z2^{*}$\\[1pt] $z1^{*}$}; &
-    &
-    \node [Optimization] (53797374656d2d4f7074696d697a6572) {System-OPT}; &
-    \node [DataInter] (4631-53797374656d2d4f7074696d697a6572) {$z2$\\[1pt] $x1^{c}$\\[1pt] $y1^{c}$\\[1pt] $y2^{c}$}; &
-    &
-    \node [DataInter] (4431-53797374656d2d4f7074696d697a6572) {$y2^{c}$}; &
-    &
-    \node [DataInter] (5f5f4a305f5f-53797374656d2d4f7074696d697a6572) {$z2$\\[1pt] $z1$\\[1pt] $y1^{c}$\\[1pt] $x1^{c}$}; &
-    &
-    \node [DataInter] (4432-53797374656d2d4f7074696d697a6572) {$y1^{c}$}; &
-    &
-    \node [DataInter] (5f5f4a315f5f-53797374656d2d4f7074696d697a6572) {$z2$\\[1pt] $z1$\\[1pt] $y2^{c}$}; \\
-    %Row 5
-    &
-    \node [DataIO] (436f6f7264696e61746f72-4631) {$f^{*}$}; &
-    &
-    \node [DataInter] (53797374656d2d4f7074696d697a6572-4631) {$f$}; &
-    \node [PostAnalysis] (4631) {F1}; &
-    &
-    &
-    &
-    &
-    &
-    &
-    &
-    \\
-    &
-    \node [DataIO] (436f6f7264696e61746f72-5375624f7074696d697a65722d30) {$x1^{*}$\\[1pt] $z1^{c,*}$\\[1pt] $z2^{c,*}$}; &
-    &
-    &
-    &
-    \node [Optimization] (5375624f7074696d697a65722d30) {SubOPT-0}; &
-    \node [DataInter] (4431-5375624f7074696d697a65722d30) {$x1$\\[1pt] $z1^{c}$\\[1pt] $z2^{c}$}; &
-    &
-    \node [DataInter] (5f5f4a305f5f-5375624f7074696d697a65722d30) {$x1$\\[1pt] $z1^{c}$\\[1pt] $z2^{c}$}; &
-    &
-    &
-    &
-    \\
-    &
-    &
-    &
-    &
-    &
-    &
-    \node [CoupledAnalysis] (4431) {D1}; &
-    \node [DataInter] (4731-4431) {$y1$}; &
-    \node [DataInter] (5f5f4a305f5f-4431) {$y1$}; &
-    &
-    &
-    &
-    \\
-    &
-    \node [DataIO] (436f6f7264696e61746f72-4731) {$g1^{*}$}; &
-    &
-    &
-    &
-    \node [DataInter] (5375624f7074696d697a65722d30-4731) {$g1$}; &
-    &
-    \node [PostAnalysis] (4731) {G1}; &
-    &
-    &
-    &
-    &
-    \\
-    &
-    \node [DataIO] (436f6f7264696e61746f72-5f5f4a305f5f) {$J0^{*,*}$}; &
-    &
-    \node [DataInter] (53797374656d2d4f7074696d697a6572-5f5f4a305f5f) {$J0^{*}$}; &
-    &
-    \node [DataInter] (5375624f7074696d697a65722d30-5f5f4a305f5f) {$J0$}; &
-    &
-    &
-    \node [PostAnalysis] (5f5f4a305f5f) {J0}; &
-    &
-    &
-    &
-    \\
-    &
-    \node [DataIO] (436f6f7264696e61746f72-5375624f7074696d697a65722d31) {$z1^{c,i2,*}$\\[1pt] $z2^{c,i2,*}$}; &
-    &
-    &
-    &
-    &
-    &
-    &
-    &
-    \node [Optimization] (5375624f7074696d697a65722d31) {SubOPT-1}; &
-    \node [DataInter] (4432-5375624f7074696d697a65722d31) {$z1^{c,i2}$\\[1pt] $z2^{c,i2}$}; &
-    &
-    \node [DataInter] (5f5f4a315f5f-5375624f7074696d697a65722d31) {$z1^{c,i2}$\\[1pt] $z2^{c,i2}$}; \\
-    %Row 11
-    &
-    &
-    &
-    &
-    &
-    &
-    &
-    &
-    &
-    &
-    \node [CoupledAnalysis] (4432) {D2}; &
-    \node [DataInter] (4732-4432) {$y2$}; &
-    \node [DataInter] (5f5f4a315f5f-4432) {$y2$}; \\
-    %Row 12
-    &
-    \node [DataIO] (436f6f7264696e61746f72-4732) {$g2^{*}$}; &
-    &
-    &
-    &
-    &
-    &
-    &
-    &
-    \node [DataInter] (5375624f7074696d697a65722d31-4732) {$g2$}; &
-    &
-    \node [PostAnalysis] (4732) {G2}; &
-    \\
-    &
-    \node [DataIO] (436f6f7264696e61746f72-5f5f4a315f5f) {$J1^{*,*}$}; &
-    &
-    \node [DataInter] (53797374656d2d4f7074696d697a6572-5f5f4a315f5f) {$J1^{*}$}; &
-    &
-    &
-    &
-    &
-    &
-    \node [DataInter] (5375624f7074696d697a65722d31-5f5f4a315f5f) {$J1$}; &
-    &
-    &
-    \node [PostAnalysis] (5f5f4a315f5f) {J1}; \\
-    %Row 14
-  };
-
-  \begin{pgfonlayer}{data}
-    \path
-    % Horizontal edges
-    (436f6f7264696e61746f72) edge [DataLine] (5375624f7074696d697a65722d31-436f6f7264696e61746f72)
-    (41) edge [DataLine] (4432-41)
-    (436f6f7264696e61746f72-53797374656d2d4f7074696d697a6572) edge [DataLine] (5f5f4a315f5f-53797374656d2d4f7074696d697a6572)
-    (436f6f7264696e61746f72-4631) edge [DataLine] (4631)
-    (436f6f7264696e61746f72-5375624f7074696d697a65722d30) edge [DataLine] (5f5f4a305f5f-5375624f7074696d697a65722d30)
-    (4431) edge [DataLine] (5f5f4a305f5f-4431)
-    (436f6f7264696e61746f72-4731) edge [DataLine] (4731)
-    (436f6f7264696e61746f72-5f5f4a305f5f) edge [DataLine] (5f5f4a305f5f)
-    (436f6f7264696e61746f72-5375624f7074696d697a65722d31) edge [DataLine] (5f5f4a315f5f-5375624f7074696d697a65722d31)
-    (4432) edge [DataLine] (5f5f4a315f5f-4432)
-    (436f6f7264696e61746f72-4732) edge [DataLine] (4732)
-    (436f6f7264696e61746f72-5f5f4a315f5f) edge [DataLine] (5f5f4a315f5f)
-    % Vertical edges
-    (436f6f7264696e61746f72) edge [DataLine] (436f6f7264696e61746f72-5f5f4a315f5f)
-    (41-436f6f7264696e61746f72) edge [DataLine] (41)
-    (53797374656d2d4f7074696d697a6572-436f6f7264696e61746f72) edge [DataLine] (53797374656d2d4f7074696d697a6572-5f5f4a315f5f)
-    (4631-53797374656d2d4f7074696d697a6572) edge [DataLine] (4631)
-    (5375624f7074696d697a65722d30-436f6f7264696e61746f72) edge [DataLine] (5375624f7074696d697a65722d30-5f5f4a305f5f)
-    (4431-41) edge [DataLine] (4431)
-    (4731-4431) edge [DataLine] (4731)
-    (5f5f4a305f5f-53797374656d2d4f7074696d697a6572) edge [DataLine] (5f5f4a305f5f)
-    (5375624f7074696d697a65722d31-436f6f7264696e61746f72) edge [DataLine] (5375624f7074696d697a65722d31-5f5f4a315f5f)
-    (4432-41) edge [DataLine] (4432)
-    (4732-4432) edge [DataLine] (4732)
-    (5f5f4a315f5f-53797374656d2d4f7074696d697a6572) edge [DataLine] (5f5f4a315f5f)
-    ;
-  \end{pgfonlayer}
-
-\end{tikzpicture}
-
-\end{document}
diff --git a/examples/scripts/sellar_problem/KDMS/Mdao_CO.kdms b/examples/scripts/sellar_problem/KDMS/Mdao_CO.kdms
new file mode 100644
index 0000000000000000000000000000000000000000..1a5c361d95aa34d67cd325732225cc55d6d2366c
Binary files /dev/null and b/examples/scripts/sellar_problem/KDMS/Mdao_CO.kdms differ
diff --git a/examples/scripts/sellar_problem/KDMS/Mdao_CO_mpg.kdms b/examples/scripts/sellar_problem/KDMS/Mdao_CO_mpg.kdms
new file mode 100644
index 0000000000000000000000000000000000000000..1d252429c6aa5987421119f60a336b244ece89af
Binary files /dev/null and b/examples/scripts/sellar_problem/KDMS/Mdao_CO_mpg.kdms differ
diff --git a/examples/scripts/sellar_problem_dev_CO.py b/examples/scripts/sellar_problem_dev_CO.py
index 969e8492cb23d872e2aed1c9f885fd20c58b7058..f781b43d3586f7516037f5db61824e33c1d0cbe7 100644
--- a/examples/scripts/sellar_problem_dev_CO.py
+++ b/examples/scripts/sellar_problem_dev_CO.py
@@ -200,10 +200,6 @@ fpg.add_function_problem_roles()
 
 # Get Mdao graphs
 mdg = fpg.get_mdg(name='mdg Sellar problem')
-
-mdg.create_dsm(file_name='Test_XDSM', include_system_vars=True, destination_folder=pdf_dir,
-               function_order=['Coordinator', 'A', 'System-Optimizer', 'F1', 'SubOptimizer-0', 'D1', 'G1', '__J0__', 'SubOptimizer-1', 'D2', 'G2', '__J1__'], keep_tex_file=True, open_pdf=True)
-
 mpg = mdg.get_mpg(name='mpg Sellar problem')
 mdg.graph['name'] = 'XDSM - ' + mdao_definition
 mdg.graph['description'] = 'Solution strategy to solve the Sellar problem using the strategy: ' \
diff --git a/examples/scripts/ssbj_update.py b/examples/scripts/ssbj_update.py
new file mode 100644
index 0000000000000000000000000000000000000000..fbc8a30c9a4e06fa0a360458695ce54b90a33099
--- /dev/null
+++ b/examples/scripts/ssbj_update.py
@@ -0,0 +1,185 @@
+# Imports
+import os
+import logging
+
+from collections import OrderedDict
+
+from kadmos.graph import FundamentalProblemGraph, load
+from kadmos.utilities.general import get_mdao_setup
+
+
+# Settings for logging
+logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.DEBUG)
+
+# List of MDAO definitions that can be wrapped around the problem
+mdao_definitions = ['MDF-GS',              # 0
+                    'MDF-J',               # 1
+                    'IDF']                 # 2
+
+# Settings for scripting
+mdao_definitions_loop_all = True      # Option for looping through all MDAO definitions
+mdao_definition_id = 2                # Option for selecting a MDAO definition (in case mdao_definitions_loop_all=False)
+
+# Settings for creating the CMDOWS files
+create_rcg_cmdows = True              # Option for creating the RCG CMDOWS file, set to False to save time
+
+# Settings for creating visualizations
+create_vis = True                     # Create visualisations
+create_rcg_vis = True                 # Create RCG visualizations, set to False after first execution to save time
+
+# Settings for loading and saving
+kb_dir = os.path.join(os.path.dirname(os.path.realpath(__file__)), '../knowledgebases')
+pdf_dir = 'ssbj/(X)DSM'
+cmdows_dir = 'ssbj/CMDOWS'
+kdms_dir = 'ssbj/KDMS'
+vistoms_dir = 'ssbj/VISTOMS'
+
+print 'Loading repository connectivity graph...'
+
+rcg = load(os.path.join(kb_dir, 'ssbj', '__cmdows__SSBJ.xml'), io_xsd_check=True)
+
+print 'Scripting RCG...'
+
+# A name and a description are added to the graph
+rcg.graph['name'] = 'RCG'
+rcg.graph['description'] = 'Repository of the super-sonic business jet test case optimization problem'
+
+# Add some (optional) organization information
+contacts = [{'attrib': {'uID': 'ivangent'}, 'name': 'Imco van Gent', 'email': 'i.vangent@tudelft.nl', 'company': 'TU Delft'},
+            {'attrib': {'uID': 'lmuller'}, 'name': 'Lukas Muller', 'email': 'l.muller@student.tudelft.nl', 'company': 'TU Delft'}]
+architects = [{'contactUID': 'ivangent'}, {'contactUID': 'lmuller'}]
+integrators = [{'contactUID': 'lmuller'}]
+rcg.graph['organization'] = OrderedDict([('contacts', contacts),
+                                         ('organigram', {'architects': architects,
+                                                         'integrators': integrators})])
+
+# Add the objective
+rcg.add_node('objective', category='function')
+rcg.add_node('/data_schema/aircraft/other/objective', category='variable', label='obj')
+rcg.add_edge('/data_schema/aircraft/other/R', 'objective')
+rcg.add_edge('objective', '/data_schema/aircraft/other/objective')
+rcg.add_equation_labels(rcg.get_function_nodes())
+rcg.add_equation('objective', '-R', 'Python')
+rcg.add_equation('objective', '-R', 'LaTeX')
+
+# Define function order for visualization (otherwise the functions will be placed randomly on the diagonal)
+functions = ['structure[main][1][1.0]',
+             'aerodynamics[main][1][1.0]',
+             'propulsion[main][1][1.0]',
+             'performance[main][1][1.0]',
+             'objective']
+
+# Create a DSM and a VISTOMS visualization of the RCG
+if create_vis and create_rcg_vis:
+    rcg.create_dsm('RCG', include_system_vars=True, function_order=functions,
+                   destination_folder=pdf_dir)
+    rcg.vistoms_create(vistoms_dir, function_order=functions)
+
+# Save CMDOWS file
+if create_rcg_cmdows:
+    rcg.save('RCG',
+             file_type='cmdows',
+             description='RCG CMDOWS file of the super-sonic business jet test case optimization problem',
+             creator='Lukas Mueller',
+             version='0.1',
+             destination_folder=cmdows_dir,
+             pretty_print=True,
+             integrity=True)
+
+
+# On to the wrapping of the MDAO architectures
+# Get iterator (all or single one)
+if not mdao_definitions_loop_all:
+    mdao_definitions = [mdao_definitions[mdao_definition_id]]
+
+for mdao_definition in mdao_definitions:
+
+    print 'Scripting ' + str(mdao_definition) + '...'
+
+    # Reset FPG
+    fpg = FundamentalProblemGraph(rcg)
+    fpg.graph['name'] = rcg.graph['name'] + ' - ' + mdao_definition + ' - FPG'
+    fpg.graph['description'] = 'Fundamental problem graph to solve the super-sonic business jet test case optimization problem using the strategy: ' \
+                               + mdao_definition + '.'
+
+    # Determine the three main settings: architecture, convergence type and unconverged coupling setting
+    mdao_architecture, convergence_type, allow_unconverged_couplings = get_mdao_setup(mdao_definition)
+
+    # Define settings of the problem formulation
+    fpg.graph['problem_formulation'] = dict()
+    fpg.graph['problem_formulation']['function_order'] = functions
+    fpg.graph['problem_formulation']['mdao_architecture'] = mdao_architecture
+    fpg.graph['problem_formulation']['convergence_type'] = convergence_type
+    fpg.graph['problem_formulation']['allow_unconverged_couplings'] = allow_unconverged_couplings
+
+    # Depending on the architecture, different additional node attributes have to be specified. This is automated here
+    # to allow direct execution of all different options.
+    if mdao_architecture in ['IDF', 'MDF']:
+        fpg.nodes['/data_schema/aircraft/other/objective']['problem_role'] = 'objective'
+        fpg.nodes['/data_schema/aircraft/geometry/tc']['problem_role'] = 'design variable'
+        fpg.nodes['/data_schema/aircraft/reference/h']['problem_role'] = 'design variable'
+        fpg.nodes['/data_schema/reference/M']['problem_role'] = 'design variable'
+        fpg.nodes['/data_schema/aircraft/geometry/AR']['problem_role'] = 'design variable'
+        fpg.nodes['/data_schema/aircraft/geometry/Lambda']['problem_role'] = 'design variable'
+        fpg.nodes['/data_schema/aircraft/geometry/Sref']['problem_role'] = 'design variable'
+        fpg.nodes['/data_schema/aircraft/geometry/lambda']['problem_role'] = 'design variable'
+        fpg.nodes['/data_schema/aircraft/geometry/section']['problem_role'] = 'design variable'
+        fpg.nodes['/data_schema/aircraft/other/Cf']['problem_role'] = 'design variable'
+        fpg.nodes['/data_schema/aircraft/other/T']['problem_role'] = 'design variable'
+        #fpg.nodes['/data_schema/aircraft/geometry/Theta']['problem_role'] = 'design variable'
+        #fpg.nodes['/data_schema/aircraft/other/L']['problem_role'] = 'design variable'
+        #fpg.nodes['/data_schema/aircraft/weight/WE']['problem_role'] = 'design variable'
+        #fpg.nodes['/data_schema/aircraft/weight/WT']['problem_role'] = 'design variable'
+        #fpg.nodes['/data_schema/reference/ESF']['problem_role'] = 'design variable'
+        #fpg.nodes['/data_schema/aircraft/other/D']['problem_role'] = 'design variable'
+        fpg.nodes['/data_schema/aircraft/other/DT']['problem_role'] = 'constraint'
+        fpg.nodes['/data_schema/aircraft/other/sigma']['problem_role'] = 'constraint'
+        fpg.nodes['/data_schema/aircraft/other/dpdx']['problem_role'] = 'constraint'
+        fpg.nodes['/data_schema/reference/Temp']['problem_role'] = 'constraint'
+
+    # Search for problem roles
+    fpg.add_function_problem_roles()
+
+    # Create a DSM visualization of the FPG
+    fpg.create_dsm(file_name='FPG_' + mdao_definition, function_order=functions, include_system_vars=True,
+                   destination_folder=pdf_dir)
+    # Create a VISTOMS visualization of the FPG (and add it to the existing directory)
+    fpg.vistoms_add(vistoms_dir, function_order=functions)
+
+    # Save the FPG as kdms
+    fpg.save('FPG_' + mdao_definition, destination_folder=kdms_dir)
+    # Save the FPG as cmdows (and do an integrity check)
+    fpg.save('FPG_' + mdao_definition, file_type='cmdows', destination_folder=cmdows_dir,
+             description='FPG CMDOWS file of the super-sonic business jet test case optimization problem',
+             creator='Imco van Gent',
+             version='0.1',
+             pretty_print=True,
+             integrity=True)
+
+    # Get Mdao graphs
+    mpg = fpg.get_mpg(name='mpg Sellar problem')
+    mdg = fpg.get_mdg(name='mdg Sellar problem')
+    mdg.graph['name'] = rcg.graph['name'] + ' - ' + mdao_definition + ' - Mdao'
+    mdg.graph['description'] = 'Solution strategy to solve the super-sonic business jet test case optimization problem using the strategy: ' \
+                               + str(mdao_architecture) + (
+                               '_' + str(convergence_type) if convergence_type else '') + '.'
+
+    # Create a DSM visualization of the Mdao
+    mdg.create_dsm(file_name='Mdao_' + mdao_definition, include_system_vars=True, destination_folder=pdf_dir,
+                   mpg=mpg)
+    # Create a VISTOMS visualization of the Mdao (and add it to the existing directory)
+    mdg.vistoms_add(vistoms_dir, mpg=mpg)
+
+    # Save the Mdao as kdms
+    mdg.save('Mdao_' + mdao_definition, destination_folder=kdms_dir, mpg=mpg)
+    # Save the Mdao as cmdows (and do an integrity check)
+    mdg.save('Mdao_' + mdao_definition, file_type='cmdows', destination_folder=cmdows_dir,
+             mpg=mpg,
+             description='Mdao CMDOWS file of the super-sonic business jet test case optimization problem',
+             creator='Imco van Gent',
+             version='0.1',
+             pretty_print=True,
+             integrity=True
+             )
+
+print 'Done!'
diff --git a/kadmos/graph/graph_data.py b/kadmos/graph/graph_data.py
index b231e657ef67f312a6731df33eedd6556ea39b81..08518f089f16e7f0b46214df0a449406097e5578 100644
--- a/kadmos/graph/graph_data.py
+++ b/kadmos/graph/graph_data.py
@@ -2529,7 +2529,7 @@ class FundamentalProblemGraph(DataGraph, KeChainMixin):
         for post_desvar in post_desvars:
             for local_function in local_functions:
                 if nx.has_path(subgraph, post_desvar, local_function):
-                    sub_level_function_dict[self.FUNCTION_ROLES[4]].append()
+                    sub_level_function_dict[self.FUNCTION_ROLES[4]].append(post_desvar)
         additional_post_couplings = []
         for post_coupling in post_couplings:
             if post_coupling not in local_functions:
@@ -2733,11 +2733,10 @@ class FundamentalProblemGraph(DataGraph, KeChainMixin):
             pass
         elif mdao_arch == graph.OPTIONS_ARCHITECTURES[8]: # CO
             coupled_functions_groups = graph.graph['problem_formulation']['coupled_functions_groups']
-            n_groups = len(coupled_functions_groups)
-            sys_opt = 'System-{}'.format(graph.OPTIMIZER_STRING)
-            sys_opt_label = 'System-{}'.format(graph.OPTIMIZER_LABEL)
-            sub_opts = ['Sub{}-{}'.format(graph.OPTIMIZER_STRING, item[0]) for item in enumerate(coupled_functions_groups)]
-            sub_opts_labels = ['Sub{}-{}'.format(graph.OPTIMIZER_LABEL, item[0]) for item in enumerate(coupled_functions_groups)]
+            sys_opt = '{}{}'.format(graph.SYS_PREFIX, graph.OPTIMIZER_STRING)
+            sys_opt_label = '{}{}'.format(graph.SYS_PREFIX, graph.OPTIMIZER_LABEL)
+            sub_opts = ['{}{}{}{}'.format(graph.SUBSYS_PREFIX, graph.OPTIMIZER_STRING, graph.SUBSYS_SUFFIX, item[0]) for item in enumerate(coupled_functions_groups)]
+            sub_opts_labels = ['{}{}{}{}'.format(graph.SUBSYS_PREFIX, graph.OPTIMIZER_LABEL, graph.SUBSYS_SUFFIX, item[0]) for item in enumerate(coupled_functions_groups)]
 
             # Determine coupling variables between coupled_function_groups (these become design variables of the system optimizer)
             group_couplings, group_couplings_groups_idx = graph.get_group_couplings()
@@ -2759,7 +2758,6 @@ class FundamentalProblemGraph(DataGraph, KeChainMixin):
             local_cnstrnt_funcs, \
             cnstrnt_vars_group_idxs, \
             cnstrnt_funcs_group_idxs = graph.determine_scope_constraint_functions(cnstrnt_vars=constraint_nodes)
-            # TODO: determine other post-coupling functions that have a dependency with the constraint function
 
             # Create dictionary of pre-desvar, post-desvar, and post-coupling functions for the system optimizer
             sys_functions_dict = graph.get_system_level_functions(global_objective_function, global_cnstrnt_funcs,
@@ -2779,11 +2777,16 @@ class FundamentalProblemGraph(DataGraph, KeChainMixin):
                 # Create dict collecting the subsystem functions dictionary
                 subsys_functions_dicts.append(subsys_functions_dict)
 
+            # TODO: Assert that function instances are not required (future functionality)
+
             # TODO: Determine the functions that require instances, add them, and adjust subsys_functions_dict accordingly
+            # sys_functions_dict, subsys_functions_dicts = graph.create_function_instances(sys_functions_dict,
+            #                                                                              subsys_functions_dicts)
 
             # Keep track of the design variables and constraints for the system level
             sys_lev_des_vars = set(global_des_vars)
             sys_lev_cnstrnts = set(global_cnstrnt_vars)
+
             # For each discipline group, localize the group, add the consistency objective function and add the
             # sub-optimizer
             for idx, subsys_functions_dict in enumerate(subsys_functions_dicts):
@@ -2813,6 +2816,7 @@ class FundamentalProblemGraph(DataGraph, KeChainMixin):
                 cof_mappings = mapping_des_vars.copy()
                 cof_mappings.update(mapping_locals)
                 group_cof_node, group_cof_obj_node = mdg.connect_consistency_objective_function(idx, cof_mappings)
+                subsys_functions_dicts[idx][mdg.FUNCTION_ROLES[2]].append(group_cof_node)
 
                 # TODO: Then (optionally) add a converger or check for the removal of feedback (NO!!)?
                 # if feedback inside the coupled group
@@ -2854,6 +2858,9 @@ class FundamentalProblemGraph(DataGraph, KeChainMixin):
             # Finally, connect the coordinator
             mdg.connect_coordinator()
 
+            # Write function_ordering to the graph
+            mdg.graph['distr_function_ordering'] = [sys_functions_dict, subsys_functions_dicts]
+
         logger.info('Composed MDG.')
 
         return mdg
@@ -3282,19 +3289,25 @@ class MdaoDataGraph(DataGraph, MdaoMixin):
 
         # Add the consistency constraint objective function (as generic function node, since it will be made as a
         # mathematical function)
-        new_function_node = '__J' + str(group_idx) + '__'
+        new_function_node = self.COF_STRING + str(group_idx) + self.COF_SUFFIX
         assert not self.has_node(new_function_node), \
             'The automatically generated function {} somehow already exists.'.format(new_function_node)
         self.add_node(new_function_node,
                       category='function',
-                      label='J' + str(group_idx),
+                      label=self.COF_LABEL + str(group_idx),
                       instance=1,
                       problem_role=self.FUNCTION_ROLES[2],  # post-coupling
                       architecture_role=self.ARCHITECTURE_ROLES_FUNS[9])  # consistency constraint function
         # Connect the variable inputs for the function
-        for var1, var2 in ccv_mappings.iteritems():
-            self.add_edge(var1, new_function_node) # TODO: Add equation_label
-            self.add_edge(var2, new_function_node) # TODO: Add equation_label
+        for idx, (var1, var2) in enumerate(ccv_mappings.iteritems()):
+            eq_lab1 = 'x{}_0'.format(idx)
+            eq_lab2 = 'x{}_1'.format(idx)
+            self.add_edge(var1, new_function_node, equation_label=eq_lab1)
+            self.add_edge(var2, new_function_node, equation_label=eq_lab2)
+            if idx == 0:
+                math_expression = '({}-{})**2'.format(eq_lab2, eq_lab1)
+            else:
+                math_expression += '+({}-{})**2'.format(eq_lab2, eq_lab1)
         # Create the output objective node of the function and connect it
         xpath_var1 = var1.split('/')
         root = xpath_var1[1]
@@ -3304,7 +3317,8 @@ class MdaoDataGraph(DataGraph, MdaoMixin):
                       label='J'+str(group_idx),
                       instance=1,
                       problem_role=self.PROBLEM_ROLES_VARS[1])  # objective
-        self.add_edge(new_function_node, new_obj_node)  # TODO: Add Python equation to the edge
+        self.add_edge(new_function_node, new_obj_node)
+        self.add_equation((new_function_node, new_obj_node), math_expression, 'Python')
         return new_function_node, new_obj_node
 
     def localize_design_variables(self, group_functions, global_des_vars, local_des_vars):
@@ -3779,7 +3793,7 @@ class MdaoDataGraph(DataGraph, MdaoMixin):
         coor = mdg.COORDINATOR_STRING
         mdao_arch = mdg.graph['problem_formulation']['mdao_architecture']
 
-        # Get the function ordering for the FPG and assign function lists accordingly.
+        # Get the monolithic function ordering from the MDG and assign function lists accordingly.
         mg_function_ordering = mdg.graph['mg_function_ordering']
         if mdg.FUNCTION_ROLES[0] in mg_function_ordering:
             pre_functions = mg_function_ordering[mdg.FUNCTION_ROLES[0]]
@@ -3861,8 +3875,22 @@ class MdaoDataGraph(DataGraph, MdaoMixin):
         elif mdao_arch == mdg.OPTIONS_ARCHITECTURES[7]:  # distributed-convergence
             pass
         elif mdao_arch == mdg.OPTIONS_ARCHITECTURES[8]:  # CO
-            pass
-
+            distr_function_ordering = mdg.graph['distr_function_ordering']
+            sys_opt = '{}{}'.format(mdg.SYS_PREFIX, mdg.OPTIMIZER_STRING)
+            sub_opts = ['{}{}{}{}'.format(mdg.SUBSYS_PREFIX, mdg.OPTIMIZER_STRING, mdg.SUBSYS_SUFFIX, item[0]) for
+                        item in enumerate(distr_function_ordering[1])]
+            sequence1 = [coor] + distr_function_ordering[0][self.FUNCTION_ROLES[3]] + [sys_opt]
+            mpg.add_process(sequence1, 0, mdg)
+            sequence2 = [sys_opt] + distr_function_ordering[0][self.FUNCTION_ROLES[2]]
+            mpg.add_process(sequence2, mpg.node[sequence1[-1]]['process_step'], mdg, end_in_iterative_node=sys_opt)
+            for idx, subgroup in enumerate(distr_function_ordering[1]):
+                sequence3 = [sys_opt, sub_opts[idx]]
+                mpg.connect_nested_iterators(sys_opt, sub_opts[idx], direction='master->slave')
+                sequence4 = [sub_opts[idx]] + subgroup[self.FUNCTION_ROLES[4]] + subgroup[self.FUNCTION_ROLES[1]] + \
+                            subgroup[self.FUNCTION_ROLES[2]]
+                mpg.add_process(sequence4, mpg.node[sequence3[-1]]['process_step'], mdg, end_in_iterative_node=sub_opts[idx])
+                mpg.connect_nested_iterators(sys_opt, sub_opts[idx])
+            mpg.connect_nested_iterators(coor, sys_opt)
 
         mpg.graph['process_hierarchy'] = mpg.get_process_hierarchy()
 
diff --git a/kadmos/graph/graph_kadmos.py b/kadmos/graph/graph_kadmos.py
index 83b86294009e3bcfe9aa8ff61623da355cfd63c8..0cde16bc5bd7d53cecfb4c2fe2a736d1f793c8ae 100644
--- a/kadmos/graph/graph_kadmos.py
+++ b/kadmos/graph/graph_kadmos.py
@@ -116,12 +116,18 @@ class KadmosGraph(nx.DiGraph, EquationMixin, VistomsMixin):
                                      'PostAnalysis',               # 8
                                      'PostAnalysis']               # 9
     CMDOWS_ARCHITECTURE_ROLE_SPLITTER = get_list_entries(ARCHITECTURE_ROLES_FUNS, 0, 1, 2, 3, 9)
+    SYS_PREFIX = 'Sys-'
+    SUBSYS_PREFIX = 'Sub-'
+    SUBSYS_SUFFIX = '-'
     COORDINATOR_STRING = 'Coordinator'
     COORDINATOR_LABEL = 'COOR'
     CONVERGER_STRING = 'Converger'
     CONVERGER_LABEL = 'CONV'
     CONSCONS_STRING = 'Gc'
     CONSCONS_LABEL = 'Gc'
+    COF_STRING = '__J'
+    COF_SUFFIX = '__'
+    COF_LABEL = 'J'
     DOE_STRING = 'DOE'
     DOE_LABEL = 'DOE'
     OPTIMIZER_STRING = 'Optimizer'
@@ -4028,7 +4034,7 @@ def _load_cmdows(file_path, io_xsd_check):
     elif rcg:
         graph = RepositoryConnectivityGraph()
     else:
-        IOError('The CMDOWS file ' + file_path + ' is missing basic elements and cannot be loaded.')
+        raise IOError('The CMDOWS file ' + file_path + ' is missing basic elements and cannot be loaded.')
 
     # Load the graph (and MPG in case this one is provided)
     mpg = graph.load_cmdows(cmdows, io_xsd_check)
diff --git a/kadmos/graph/graph_process.py b/kadmos/graph/graph_process.py
index 1ab689281fc0c003cb925d391178511636ebe52d..7c6183d8d6a4cc67247108fb11e7ebddc397fe53 100644
--- a/kadmos/graph/graph_process.py
+++ b/kadmos/graph/graph_process.py
@@ -209,43 +209,145 @@ class MdaoProcessGraph(ProcessGraph):
         :rtype: None
         """
 
-        # Get function ordering of MDAO graph and establish diagonal order list
-        mg_function_ordering = self.graph['mg_function_ordering']
-        diagonal_order = self.find_all_nodes(attr_cond=['architecture_role', '==',
-                                                        self.ARCHITECTURE_ROLES_FUNS[0]])  # coordinator
-
-        # Append pre-coupling functions
-        if self.FUNCTION_ROLES[0] in mg_function_ordering:
-            diagonal_order.extend(mg_function_ordering[self.FUNCTION_ROLES[0]])
-
-        # Append pre-desvars functions
-        if self.FUNCTION_ROLES[3] in mg_function_ordering:
-            diagonal_order.extend(mg_function_ordering[self.FUNCTION_ROLES[3]])
-
-        # Append optimizer or DOE block
-        diagonal_order.extend(self.find_all_nodes(attr_cond=['architecture_role', '==',
-                                                             self.ARCHITECTURE_ROLES_FUNS[1]]))  # optimizer
-        diagonal_order.extend(self.find_all_nodes(attr_cond=['architecture_role', '==',
-                                                             self.ARCHITECTURE_ROLES_FUNS[3]]))  # doe
-
-        # Append post-desvars functions
-        if self.FUNCTION_ROLES[4] in mg_function_ordering:
-            diagonal_order.extend(mg_function_ordering[self.FUNCTION_ROLES[4]])
-
-        # Append converger block
-        diagonal_order.extend(self.find_all_nodes(attr_cond=['architecture_role', '==',
-                                                             self.ARCHITECTURE_ROLES_FUNS[2]]))  # converger
-
-        # Append coupled functions
-        if self.FUNCTION_ROLES[1] in mg_function_ordering:
-            diagonal_order.extend(mg_function_ordering[self.FUNCTION_ROLES[1]])
-
-        # Append post-coupling functions
-        if self.FUNCTION_ROLES[2] in mg_function_ordering:
-            diagonal_order.extend(mg_function_ordering[self.FUNCTION_ROLES[2]])
+        # TODO: Update this function to only one function_ordering for both monolithic and distributed architectures
+
+        if 'distr_function_ordering' not in self.graph:
+            # Get function ordering of MDAO graph and establish diagonal order list
+            mg_function_ordering = self.graph['mg_function_ordering']
+            diagonal_order = self.find_all_nodes(attr_cond=['architecture_role', '==',
+                                                            self.ARCHITECTURE_ROLES_FUNS[0]])  # coordinator
+
+            # Append pre-coupling functions
+            if self.FUNCTION_ROLES[0] in mg_function_ordering:
+                diagonal_order.extend(mg_function_ordering[self.FUNCTION_ROLES[0]])
+
+            # Append pre-desvars functions
+            if self.FUNCTION_ROLES[3] in mg_function_ordering:
+                diagonal_order.extend(mg_function_ordering[self.FUNCTION_ROLES[3]])
+
+            # Append optimizer or DOE block
+            diagonal_order.extend(self.find_all_nodes(attr_cond=['architecture_role', '==',
+                                                                 self.ARCHITECTURE_ROLES_FUNS[1]]))  # optimizer
+            diagonal_order.extend(self.find_all_nodes(attr_cond=['architecture_role', '==',
+                                                                 self.ARCHITECTURE_ROLES_FUNS[3]]))  # doe
+
+            # Append post-desvars functions
+            if self.FUNCTION_ROLES[4] in mg_function_ordering:
+                diagonal_order.extend(mg_function_ordering[self.FUNCTION_ROLES[4]])
+
+            # Append converger block
+            diagonal_order.extend(self.find_all_nodes(attr_cond=['architecture_role', '==',
+                                                                 self.ARCHITECTURE_ROLES_FUNS[2]]))  # converger
+
+            # Append coupled functions
+            if self.FUNCTION_ROLES[1] in mg_function_ordering:
+                diagonal_order.extend(mg_function_ordering[self.FUNCTION_ROLES[1]])
+
+            # Append post-coupling functions
+            if self.FUNCTION_ROLES[2] in mg_function_ordering:
+                diagonal_order.extend(mg_function_ordering[self.FUNCTION_ROLES[2]])
+
+            for diag_pos, node in enumerate(diagonal_order):
+                self.nodes[node]['diagonal_position'] = diag_pos
+        else:
+            mg_function_ordering = self.graph['distr_function_ordering']
+            syslevel_ordering = mg_function_ordering[0]
+            subsyslevel_orderings = mg_function_ordering[1]
+            diagonal_order = self.find_all_nodes(attr_cond=['architecture_role', '==',
+                                                            self.ARCHITECTURE_ROLES_FUNS[0]])  # coordinator
+
+            # Append system-level pre-coupling functions
+            if self.FUNCTION_ROLES[0] in syslevel_ordering:
+                diagonal_order.extend(syslevel_ordering[self.FUNCTION_ROLES[0]])
+
+            # Append system-level pre-desvars functions
+            if self.FUNCTION_ROLES[3] in syslevel_ordering:
+                diagonal_order.extend(syslevel_ordering[self.FUNCTION_ROLES[3]])
+
+            # Append system level optimizer and/or DOE block
+            opts = self.find_all_nodes(attr_cond=['architecture_role', '==', self.ARCHITECTURE_ROLES_FUNS[1]]) # optimizer
+            if len(opts) > 1:
+                sys_opt = [item for item in opts if self.SYS_PREFIX in item]
+                assert len(sys_opt) == 1, '{} system optimizers found, one expected.'.format(len(sys_opt))
+                opts = sys_opt
+            diagonal_order.extend(opts)
+            does = self.find_all_nodes(attr_cond=['architecture_role', '==', self.ARCHITECTURE_ROLES_FUNS[3]])  # doe
+            if len(does) > 1:
+                sys_doe = [item for item in does if self.SYS_PREFIX in item]
+                assert len(sys_doe) == 1, '{} system DOE(s) found, one expected.'.format(len(sys_doe))
+                does = sys_doe
+            diagonal_order.extend(does)
+
+            # Append system-level post-desvars functions
+            if self.FUNCTION_ROLES[4] in syslevel_ordering:
+                diagonal_order.extend(syslevel_ordering[self.FUNCTION_ROLES[4]])
+
+            # Append system-level converger block
+            convs = self.find_all_nodes(attr_cond=['architecture_role', '==', self.ARCHITECTURE_ROLES_FUNS[2]])  # converger
+            if len(convs) > 1:
+                sys_conv = [item for item in convs if self.SYS_PREFIX in item]
+                assert len(sys_conv) == 1, '{} system convergers found, one expected.'.format(len(sys_conv))
+                convs = sys_conv
+            diagonal_order.extend(convs)  # converger
+
+            # Append system-level coupled functions
+            if self.FUNCTION_ROLES[1] in syslevel_ordering:
+                diagonal_order.extend(syslevel_ordering[self.FUNCTION_ROLES[1]])
+
+            # Append system-level post-coupling functions
+            if self.FUNCTION_ROLES[2] in syslevel_ordering:
+                diagonal_order.extend(syslevel_ordering[self.FUNCTION_ROLES[2]])
+
+            # Append sublevel functions here
+            for idx, subsyslevel_ord in enumerate(subsyslevel_orderings):
+                # Append subsystem-level pre-coupling functions
+                if self.FUNCTION_ROLES[0] in subsyslevel_ord:
+                    diagonal_order.extend(subsyslevel_ord[self.FUNCTION_ROLES[0]])
+
+                # Append subsystem-level pre-desvars functions
+                if self.FUNCTION_ROLES[3] in subsyslevel_ord:
+                    diagonal_order.extend(subsyslevel_ord[self.FUNCTION_ROLES[3]])
+
+                # Append subsytem-level optimizer and/or DOE block
+                opts = self.find_all_nodes(
+                    attr_cond=['architecture_role', '==', self.ARCHITECTURE_ROLES_FUNS[1]])  # optimizer
+                if len(opts) > 1:
+                    sys_opt = [item for item in opts if self.SUBSYS_SUFFIX in item and self.SUBSYS_SUFFIX+str(idx) in item]
+                    assert len(sys_opt) == 1, '{} subsystem optimizers found, one expected.'.format(len(sys_opt))
+                    opts = sys_opt
+                diagonal_order.extend(opts)
+                does = self.find_all_nodes(
+                    attr_cond=['architecture_role', '==', self.ARCHITECTURE_ROLES_FUNS[3]])  # doe
+                if len(does) > 1:
+                    sys_doe = [item for item in does if self.SUBSYS_SUFFIX in item and self.SUBSYS_SUFFIX+str(idx) in item]
+                    assert len(sys_doe) == 1, '{} subsystem DOEs found, one expected.'.format(len(sys_doe))
+                    does = sys_doe
+                diagonal_order.extend(does)
+
+                # Append subsytem-level post-desvars functions
+                if self.FUNCTION_ROLES[4] in subsyslevel_ord:
+                    diagonal_order.extend(subsyslevel_ord[self.FUNCTION_ROLES[4]])
+
+                # Append subsytem-level converger block
+                convs = self.find_all_nodes(
+                    attr_cond=['architecture_role', '==', self.ARCHITECTURE_ROLES_FUNS[2]])  # converger
+                if len(convs) > 1:
+                    sys_conv = [item for item in convs if self.SUBSYS_SUFFIX in item and self.SUBSYS_SUFFIX+str(idx) in item]
+                    assert len(sys_conv) == 1, '{} subsystem convergers found, one expected.'.format(len(sys_conv))
+                    convs = sys_conv
+                diagonal_order.extend(convs)  # converger
+
+                # Append subsytem-level coupled functions
+                if self.FUNCTION_ROLES[1] in subsyslevel_ord:
+                    diagonal_order.extend(subsyslevel_ord[self.FUNCTION_ROLES[1]])
+
+                # Append subsytem-level post-coupling functions
+                if self.FUNCTION_ROLES[2] in subsyslevel_ord:
+                    diagonal_order.extend(subsyslevel_ord[self.FUNCTION_ROLES[2]])
+
+            for diag_pos, node in enumerate(diagonal_order):
+                self.nodes[node]['diagonal_position'] = diag_pos
 
-        for diag_pos, node in enumerate(diagonal_order):
-            self.nodes[node]['diagonal_position'] = diag_pos
 
         return
 
@@ -339,14 +441,18 @@ class MdaoProcessGraph(ProcessGraph):
                 # End in iterative node or simply end function
                 if end_in_iterative_node:
                     [self.add_edge(node_1, end_in_iterative_node, process_step=process_step) for node_1 in nodes_1]
-                    self.nodes[end_in_iterative_node]['converger_step'] = process_step
+                    if 'converger_step' not in self.nodes[end_in_iterative_node]:
+                        self.nodes[end_in_iterative_node]['converger_step'] = process_step
+                    else:
+                        if process_step > self.nodes[end_in_iterative_node]['converger_step']:
+                            self.nodes[end_in_iterative_node]['converger_step'] = process_step
             # Increment process step
             process_step += 1
             coupling_matrix_1 = coupling_matrix_2
 
         return
 
-    def connect_nested_iterators(self, master, slave):
+    def connect_nested_iterators(self, master, slave, direction='slave->master'):
         """Method to connect a slave iterator to a master iterator in a nested configuration.
 
         An example is if a converger inside an optimizer in MDF needs to be linked back.
@@ -356,13 +462,22 @@ class MdaoProcessGraph(ProcessGraph):
         :param slave: lower iterator node in the nested configuration
         :type slave: basestring
         """
-
-        assert self.has_node(master), 'Node %s not present in the graph.' % master
-        assert self.has_node(slave), 'Node %s not present in the graph.' % slave
-        assert 'converger_step' in self.nodes[slave], 'Slave node %s needs to have a converger_step.' % slave
-        self.add_edge(slave, master, process_step=self.nodes[slave]['converger_step'] + 1)
-        self.nodes[master]['converger_step'] = self.nodes[slave]['converger_step'] + 1
-
+        dir_options = ['slave->master', 'master->slave']
+        assert direction in dir_options, 'direction options are {} and {}.'.format(dir_options[0], dir_options[1])
+        assert self.has_node(master), 'Node {} not present in the graph.'.format(master)
+        assert self.has_node(slave), 'Node {} not present in the graph.'.format(slave)
+        if direction == 'slave->master':
+            assert 'converger_step' in self.nodes[slave], 'Slave node %s needs to have a converger_step.' % slave
+            self.add_edge(slave, master, process_step=self.nodes[slave]['converger_step'] + 1)
+            if 'converger_step' not in self.nodes[master]:
+                self.nodes[master]['converger_step'] = self.nodes[slave]['converger_step'] + 1
+            else:
+                if self.nodes[slave]['converger_step'] + 1 > self.nodes[master]['converger_step']:
+                    self.nodes[master]['converger_step'] = self.nodes[slave]['converger_step'] + 1
+        else:
+            assert 'process_step' in self.nodes[master], 'Master node {} needs to have a process_step.'.format(master)
+            self.add_edge(master, slave, process_step=self.nodes[master]['process_step'] + 1)
+            self.nodes[slave]['process_step'] = self.nodes[master]['process_step'] + 1
         return
 
     def get_node_text(self, node):
diff --git a/kadmos/graph/mixin_equation.py b/kadmos/graph/mixin_equation.py
index 76d87422e1aec9a07ade64ac6c00c9236ab57ea9..81b284fb8e75ac20e1acbb71b90b84f9431c66d8 100644
--- a/kadmos/graph/mixin_equation.py
+++ b/kadmos/graph/mixin_equation.py
@@ -102,9 +102,6 @@ class EquationMixin(object):
         else:
             raise IOError('Invalid setting label_method.')
 
-
-
-
         # Make the label valid
         for char in self._get_equation_chars(language=language):
             label = label.replace(char, '')