Variable complexity problem

Former-commit-id: 96bed16e959f6abf19d6e8c97aff124e104b8cdb

examples/scripts/variable_complexity_problem.py

+# Imports
+import logging
+import numpy as np
+import random
+
+import networkx as nx
+
+from collections import OrderedDict
+
+from kadmos.graph import RepositoryConnectivityGraph, FundamentalProblemGraph
+from kadmos.utilities.general import get_mdao_setup
+
+
+# Settings for logging
+logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.WARNING)
+
+# List of MDAO definitions that can be wrapped around the problem
+mdao_definitions = ['unconverged-MDA-J',    # 0
+                    'unconverged-MDA-GS',   # 1
+                    'unconverged-DOE-GS',   # 2
+                    'unconverged-DOE-J',    # 3
+                    'converged-DOE-GS',     # 4
+                    'converged-DOE-J',      # 5
+                    'converged-MDA-J',      # 6
+                    'converged-MDA-GS',     # 7
+                    'MDF-GS',               # 8
+                    'MDF-J',                # 9
+                    'IDF']                  # 10
+
+# Settings for scripting
+mdao_definitions_loop_all = False       # Option for looping through all MDAO definitions
+mdao_definition_id = 8                 # Option for selecting a MDAO definition (in case mdao_definitions_loop_all=False)
+
+# Settings for saving
+pdf_dir = 'variable_complexity_problem/(X)DSM'
+cmdows_dir = 'variable_complexity_problem/CMDOWS'
+kdms_dir = 'variable_complexity_problem/KDMS'
+vistoms_dir = 'variable_complexity_problem/VISTOMS'
+
+# Settings for variable complexity problem
+n_disciplines = 3
+n_global_var = 2
+n_local_var = [3, 0, 1]
+n_coupling_var = [3, 1, 2]
+n_global_constraints = 3
+
+B = np.array([[1,  0,  0,  0,  3,  4],
+              [0,  5,  0,  0,  7,  8],
+              [0,  0,  9,  0,  11, 12],
+              [13, 14, 15, 16, 17, 18],
+              [0,  0,  0,  22, 23, 0],
+              [0,  0,  0,  27, 0,  28]])
+C = np.array([[0,  30],
+              [0,  32],
+              [0,  0],
+              [35, 36],
+              [0,  0],
+              [0,  0]])
+D = np.array([[41, 42, 43, 0],
+              [44, 45, 46, 0],
+              [47, 48, 49, 0],
+              [0,  0,  0,  0],
+              [0,  0,  0,  50],
+              [0,  0,  0,  51]])
+
+H = np.array([[1, 0],
+              [0, 2],
+              [3, 0]])
+I = np.array([[0, 4, 0, 0],
+              [0, 0, 5, 0],
+              [7, 0, 0, 6]])
+J = np.array([[4, 0, 0, 0, 0, 9],
+              [0, 5, 0, 0, 8, 0],
+              [0, 0, 6, 7, 0, 0]])
+s = np.array([1, 2, 3])
+
+print 'Scripting RCG...'
+
+# A new repository connectivity graph is defined to describe the general problem
+rcg = RepositoryConnectivityGraph(name='Variable complexity problem graph')
+# A description is added to the graph
+rcg.graph['description'] = 'Repository graph of the variable complexity problem'
+
+# All function nodes are defined
+for discipline in range(n_disciplines):
+    rcg.add_node('D{0}'.format(discipline+1),
+                 category='function')
+rcg.add_node('F',
+             category='function')
+for constraint in range(n_global_constraints):
+    rcg.add_node('G{0}'.format(constraint+1),
+                 category='function')
+
+# Global variable nodes are defined
+for global_var in range(n_global_var):
+    rcg.add_node('/data_schema/x0{0}'.format(global_var+1),
+                 category='variable',
+                 label='x0{0}'.format(global_var+1))
+# Local variable nodes are defined
+for discipline in range(n_disciplines):
+    for local_var in range(n_local_var[discipline]):
+        rcg.add_node('/data_schema/x{0}{1}'.format(discipline+1, local_var+1),
+                     category='variable',
+                     label='x{0}{1}'.format(discipline+1, local_var+1))
+
+# Coupling variable nodes are defined
+for discipline in range(n_disciplines):
+    for coupling_var in range(n_coupling_var[discipline]):
+        rcg.add_node('/data_schema/y{0}{1}'.format(discipline+1, coupling_var+1),
+                     category='variable',
+                     label='y{0}{1}'.format(discipline+1, coupling_var+1))
+
+# Objective variable node is defined
+rcg.add_node('/data_schema/f',
+             category='variable',
+             label='f')
+
+# Constraint variable nodes are defined
+for constraint in range(n_global_constraints):
+    rcg.add_node('/data_schema/g{0}'.format(constraint+1),
+                 category='variable',
+                 label='g{0}'.format(constraint+1))
+
+# Edges between global variables and function nodes are defined
+for global_var in range(n_global_var):
+    for discipline in range(n_disciplines):
+        values = C[sum(n_coupling_var[:discipline]):sum(n_coupling_var[:discipline+1]), global_var]
+        if not all(v == 0 for v in values):
+            rcg.add_edge('/data_schema/x0{0}'.format(global_var+1), 'D{0}'.format(discipline+1))
+    rcg.add_edge('/data_schema/x0{0}'.format(global_var+1), 'F')
+    for constraint in range(n_global_constraints):
+        # if H[constraint][global_var] != 0:
+        rcg.add_edge('/data_schema/x0{0}'.format(global_var+1), 'G{0}'.format(constraint+1))
+
+# Edges between local variables and function nodes are defined
+for discipline in range(n_disciplines):
+    for local_var in range(n_local_var[discipline]):
+        rcg.add_edge('/data_schema/x{0}{1}'.format(discipline+1, local_var+1), 'D{0}'.format(discipline+1))
+        rcg.add_edge('/data_schema/x{0}{1}'.format(discipline+1, local_var+1), 'F')
+        for constraint in range(n_global_constraints):
+            # if I[constraint][sum(n_local_var[:discipline])+local_var] != 0:
+            rcg.add_edge('/data_schema/x{0}{1}'.format(discipline+1, local_var+1), 'G{0}'.format(constraint+1))
+
+# Edges between coupling variables and function nodes are defined
+for coupling_var_disc in range(n_disciplines):
+    for coupling_var in range(n_coupling_var[coupling_var_disc]):
+        for discipline in range(n_disciplines):
+            values = B[sum(n_coupling_var[:discipline]):sum(n_coupling_var[:discipline + 1]), sum(n_coupling_var[:coupling_var_disc])+coupling_var]
+            if not discipline == coupling_var_disc and not all(v == 0 for v in values):
+                rcg.add_edge('/data_schema/y{0}{1}'.format(coupling_var_disc+1, coupling_var+1), 'D{0}'.format(discipline+1))
+        rcg.add_edge('/data_schema/y{0}{1}'.format(coupling_var_disc+1, coupling_var+1), 'F')
+        for constraint in range(n_global_constraints):
+            if J[constraint][sum(n_coupling_var[:coupling_var_disc])+coupling_var] != 0:
+                rcg.add_edge('/data_schema/y{0}{1}'.format(coupling_var_disc+1, coupling_var+1), 'G{0}'.format(constraint+1))
+
+# Edges between function nodes and coupling variables are defined
+for discipline in range(n_disciplines):
+    for coupling_var in range(n_coupling_var[discipline]):
+        rcg.add_edge('D{0}'.format(discipline+1), '/data_schema/y{0}{1}'.format(discipline+1, coupling_var+1))
+rcg.add_edge('F', '/data_schema/f')
+for constraint in range(n_global_constraints):
+    rcg.add_edge('G{0}'.format(constraint+1), '/data_schema/g{0}'.format(constraint+1))
+
+# Add equations
+rcg.add_equation_labels(rcg.get_function_nodes(), labeling_method='node_id')
+
+# Add discipline analysis equations
+for discipline in range(n_disciplines):
+    for output_var in range(n_coupling_var[discipline]):
+        equation = ""
+        for global_var in range(n_global_var):
+            if C[sum(n_coupling_var[:discipline])+output_var][global_var] != 0:
+                equation += '-{0}*x0{1}'.format(C[sum(n_coupling_var[:discipline])+output_var][global_var], global_var+1)
+        for local_var_disc in range(n_disciplines):
+            for local_var in range(n_local_var[local_var_disc]):
+                if D[sum(n_coupling_var[:discipline])+output_var][sum(n_local_var[:local_var_disc])+local_var] != 0:
+                    equation += '-{0}*x{1}{2}'.format(D[sum(n_coupling_var[:discipline])+output_var][sum(n_local_var[:local_var_disc])+local_var], local_var_disc+1, local_var+1)
+        for coupling_var_disc in range(n_disciplines):
+            for coupling_var in range(n_coupling_var[coupling_var_disc]):
+                if B[sum(n_coupling_var[:discipline])+output_var][sum(n_coupling_var[:coupling_var_disc])+coupling_var] != 0 and (discipline, output_var) != (coupling_var_disc, coupling_var):
+                    equation += '-{0}*y{1}{2}'.format(B[sum(n_coupling_var[:discipline])+output_var][sum(n_coupling_var[:coupling_var_disc])+coupling_var], coupling_var_disc+1, coupling_var+1)
+        if B[sum(n_coupling_var[:discipline])+output_var][sum(n_coupling_var[:discipline])+output_var] != 1:
+            equation = '({0})/{1}.'.format(equation, B[sum(n_coupling_var[:discipline])+output_var][sum(n_coupling_var[:discipline])+output_var])
+        rcg.add_equation(['D{0}'.format(discipline+1), '/data_schema/y{0}{1}'.format(discipline+1, output_var+1)], equation, 'Python')
+        rcg.add_equation(['D{0}'.format(discipline+1), '/data_schema/y{0}{1}'.format(discipline+1, output_var+1)], equation, 'LaTeX')
+        print equation
+
+# Add objective function equation
+objective = ""
+for global_var in range(n_global_var):
+    objective += '+x0{0}'.format(global_var+1)
+for discipline in range(n_disciplines):
+    for local_var in range(n_local_var[discipline]):
+        objective += '+x{0}{1}'.format(discipline+1, local_var+1)
+for discipline in range(n_disciplines):
+    for coupling_var in range(n_coupling_var[discipline]):
+        objective += '+y{0}{1}'.format(discipline+1, coupling_var+1)
+rcg.add_equation('F', '({0})**3'.format(objective), 'Python')
+rcg.add_equation('F', '({0})^3'.format(objective), 'LaTeX')
+print 'objective = ', '({0})^3'.format(objective)
+
+# Add constraint function equations
+for constraint in range(n_global_constraints):
+    constraint_eq = ""
+    for global_var in range(n_global_var):
+        constraint_eq += '+x0{0}*x0{0}'.format(global_var+1)
+        if H[constraint][global_var] != 0:
+            constraint_eq += '+{0}*x0{1}'.format(H[constraint][global_var], global_var+1)
+    for discipline in range(n_disciplines):
+        for local_var in range(n_local_var[discipline]):
+            constraint_eq += '+x{0}{1}*x{0}{1}'.format(discipline+1, local_var+1)
+            if I[constraint][sum(n_local_var[:discipline])+local_var] != 0:
+                constraint_eq += '+{0}*x{1}{2}'.format(I[constraint][sum(n_local_var[:discipline])+local_var], discipline+1, local_var+1)
+    for discipline in range(n_disciplines):
+        for coupling_var in range(n_coupling_var[discipline]):
+            if J[constraint][sum(n_coupling_var[:discipline])+coupling_var] != 0:
+                constraint_eq += '+{0}*y{1}{2}'.format(J[constraint][sum(n_coupling_var[:discipline])+coupling_var], discipline+1, coupling_var+1)
+    constraint_eq += '-{0}'.format(s[constraint])
+    print constraint_eq
+    rcg.add_equation('G{0}'.format(constraint+1), constraint_eq, 'Python')
+    rcg.add_equation('G{0}'.format(constraint+1), constraint_eq, 'LaTeX')
+
+function_order = []
+for discipline in range(n_disciplines):
+    function_order += ['D{0}'.format(discipline+1)]
+function_order += ['F']
+for constraint in range(n_global_constraints):
+    function_order += ['G{0}'.format(constraint+1)]
+
+# Create a DSM visualization of the RCG
+rcg.create_dsm(file_name='RCG', function_order=function_order, include_system_vars=True, destination_folder=pdf_dir)
+# Create a VISTOMS visualization of the RCG
+rcg.vistoms_create(vistoms_dir, function_order=function_order, compress=False)
+
+# Save the RCG as kdms
+rcg.save('RCG', destination_folder=kdms_dir)
+# Save the RCG as cmdows (and do an integrity check)
+rcg.save('RCG', file_type='cmdows', destination_folder=cmdows_dir,
+         description='RCG CMDOWS file of the variable complexity problem',
+         creator='Anne-Liza Bruggeman',
+         version='0.1',
+         pretty_print=True,
+         integrity=True)
+
+print 'Scripting initial FPG...'
+
+# An initial fundamental problem graph is created based on the rcg
+fpg_initial = FundamentalProblemGraph(rcg)
+
+# Determine function order
+new_graph = rcg.cleancopy()
+new_graph = new_graph.get_function_graph()
+cycles = list(nx.simple_cycles(new_graph))
+coupled_nodes = list(set(function_id for functions in cycles for function_id in functions))
+new_graph.add_node('super_node', category='function')
+for function_id in coupled_nodes:
+    new_graph = nx.contracted_nodes(new_graph, 'super_node', function_id)
+new_graph.remove_edges_from(new_graph.selfloop_edges())
+function_order = nx.topological_sort(new_graph)
+function_order[function_order.index('super_node'):function_order.index('super_node')+1] = coupled_nodes
+print function_order
+
+
+function_order = ['D1', 'D2', 'D3', 'F', 'G1', 'G2', 'G3']
+#function_order = ['D1', 'D11', 'D4', 'D3', 'D2', 'D9', 'D10', 'D6', 'D7', 'D5', 'D8', 'F', 'G1', 'G2', 'G3']
+
+# 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) + '...'
+
+    # Determine the three main settings: architecture, convergence type and unconverged coupling setting
+    mdao_architecture, convergence_type, allow_unconverged_couplings = get_mdao_setup(mdao_definition)
+
+    # Reset FPG
+    fpg = fpg_initial.cleancopy()
+    fpg.graph['name'] = rcg.graph['name'] + ' - ' + mdao_definition + ' - FPG'
+    fpg.graph['description'] = 'Fundamental problem graph for solving the variable complexity problem using the strategy: ' \
+                               + mdao_definition + '.'
+
+    # Define settings of the problem formulation
+    fpg.graph['problem_formulation'] = dict()
+    fpg.graph['problem_formulation']['function_order'] = function_order
+    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
+    if mdao_architecture in ['unconverged-DOE', 'converged-DOE']:
+        fpg.graph['problem_formulation']['doe_settings'] = dict()
+        fpg.graph['problem_formulation']['doe_settings']['doe_method'] = 'Custom design table'
+        if fpg.graph['problem_formulation']['doe_settings']['doe_method'] in ['Latin hypercube design',
+                                                                              'Monte Carlo design']:
+            fpg.graph['problem_formulation']['doe_settings']['doe_seed'] = 6
+            fpg.graph['problem_formulation']['doe_settings']['doe_runs'] = 5
+        elif fpg.graph['problem_formulation']['doe_settings']['doe_method'] in ['Full factorial design']:
+            fpg.graph['problem_formulation']['doe_settings']['doe_runs'] = 5
+
+    # 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.mark_as_objective('/data_schema/f')
+        for constraint in range(n_global_constraints):
+            fpg.mark_as_constraint('/data_schema/g{0}'.format(constraint+1), '<=', 0.0)
+        for global_var in range(n_global_var):
+            fpg.mark_as_design_variable('/data_schema/x0{0}'.format(global_var+1))
+        for discipline in range(n_disciplines):
+            for local_var in range(n_local_var[discipline]):
+                fpg.mark_as_design_variable('/data_schema/x{0}{1}'.format(discipline+1, local_var+1))
+    elif mdao_architecture in ['unconverged-MDA', 'converged-MDA']:
+        fpg.mark_as_qois(['/data_schema/f'])
+        for discipline in range(n_disciplines):
+            for coupling_var in range(n_coupling_var[discipline]):
+                fpg.mark_as_qois(['/data_schema/y{0}{1}'.format(discipline+1, coupling_var+1)])
+        for constraint in range(n_global_constraints):
+            fpg.mark_as_qois(['/data_schema/g{0}'.format(constraint+1)])
+    elif mdao_architecture in ['unconverged-DOE', 'converged-DOE']:
+        fpg.mark_as_qois(['/data_schema/f'])
+        for constraint in range(n_global_constraints):
+            fpg.mark_as_qois(['/data_schema/g{0}'.format(constraint+1)])
+        if fpg.graph['problem_formulation']['doe_settings']['doe_method'] == 'Custom design table':
+            for global_var in range(n_global_var):
+                fpg.mark_as_design_variable('/data_schema/x0{0}'.format(global_var+1), samples=[0.0, 0.1, 0.5, 0.75])
+            for discipline in range(n_disciplines):
+                for local_var in range(n_local_var[discipline]):
+                    fpg.mark_as_design_variable('/data_schema/x{0}{1}'.format(discipline+1, local_var+1), samples=[0.0, 0.1, 0.5, 0.75])
+        else:
+            for global_var in range(n_global_var):
+                fpg.mark_as_design_variable('/data_schema/x0{0}'.format(global_var+1), lower_bound=-10., upper_bound=10)
+            for discipline in range(n_disciplines):
+                for local_var in range(n_local_var[discipline]):
+                    fpg.mark_as_design_variable('/data_schema/x{0}{1}'.format(discipline+1, local_var+1), lower_bound=-10., upper_bound=10)
+    #if mdao_architecture == 'IDF':
+        #fpg.node['/data_schema/analyses/y1']['valid_ranges'] = {'limit_range': {'minimum': -100., 'maximum': 100.}}
+        #fpg.node['/data_schema/analyses/y2']['valid_ranges'] = {'limit_range': {'minimum': -100., 'maximum': 100.}}
+
+    # 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=function_order, 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=function_order)
+
+    # 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 variable complexity problem',
+             creator='Anne-Liza Bruggeman',
+             version='0.1',
+             pretty_print=True,
+             integrity=True)
+
+    # Get Mdao graphs
+    mpg = fpg.get_mpg(name='mpg variable complexity problem')
+    mdg = fpg.get_mdg(name='mdg variable complexity problem')
+    mdg.graph['name'] = rcg.graph['name'] + ' - ' + mdao_definition + ' - Mdao'
+    mdg.graph['description'] = 'Solution strategy to solve the variable complexity problem using the strategy: ' \
+                               + str(mdao_architecture) + ('_' + str(convergence_type) if convergence_type else '') + '.'
+
+    # Copy mathematical functions to additional output edges (for MDF and IDF additional output edges are created)
+    # TODO: This should be automated in the KADMOS method get_mdg()
+    if mdao_architecture in ['IDF', 'MDF']:
+        mdg.edge['F']['/data_schema/architectureNodes/finalOutputVariables/data_schemaCopy/f']['equations'] = mdg.edge['F']['/data_schema/f']['equations']
+        for constraint in range(n_global_constraints):
+            mdg.edge['G{0}'.format(constraint+1)]['/data_schema/architectureNodes/finalOutputVariables/data_schemaCopy/g{0}'.format(constraint+1)]['equations'] = mdg.edge['G{0}'.format(constraint+1)]['/data_schema/g{0}'.format(constraint+1)]['equations']
+
+    # 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)
+
+    # TODO: This should be automated
+    mdg.add_equation_labels(mdg.get_function_nodes(), labeling_method='node_id')
+
+    mdg.save('Mdao_'+mdao_definition, file_type='cmdows', destination_folder=cmdows_dir,
+             mpg=mpg,
+             description='Mdao CMDOWS file of the variable complexity problem',
+             creator='Anne-Liza Bruggeman',
+             version='0.1',
+             pretty_print=True,
+             integrity=True
+             )
+
+for discipline in range(n_disciplines):
+    if mpg.node['D{0}'.format(discipline+1)]['problem_role'] == 'coupled':
+        print 'D{0}'.format(discipline+1)
+
+print 'Done!'
\ No newline at end of file