I am attempting to implement ATC architecture into OpenMDAO v0.13.0. I've posted before about bugs and snags I've run into and have received help, but I'm still having issues. I cannot find any way to control iteration loops and the order in which sub problems are solved.
I am trying to implement this MATLAB code into equivalent OpenMDAO:

% begin ATC iterations
while conv == 0
      % solve P0
      [T11,T12,x0] = P0(x_old,R11,R12,v,w);

      % solve P11
      [R11,y11,x11] = P11(x_old,T11,y12,v,w);

      % solve P12
      [R12,y12,x12] = P12(x_old,T12,y11,v,w);

      % check convergence 
      % new vector of optimization variables
      x_new = [T11 T12 x0 R11 y11 x11 R12 y12 x12];

      % consistency constraints
      c = [T11-R11 T12-R12 y12-y11]
      c_norm(k+1) = norm(c);

      % choose criterion
      if c_norm(k+1) < ATC_eps
      % if (norm(x_new - x_old))/norm(x_old) < ATC_eps
         conv = 1;
         disp('convergence achieved')
      end
      if k > 1000
         conv = 1;
         disp('maximum number of ATC iterations reached')
      end

      %prepare for reiteration
      x_old = x_new;

      %update penalty weights
      v = v + 2*w.*w.*c;
      w = beta*w;

      k = k + 1;
end

Where P0(), P11(), and P12() are functions applying fconmin to the sub problems with the associated constraints.
My code seems to be doing this process in reverse order, and is giving incorrect results. Any help you could provide on implementing distributed architectures, architectures with penalty weights, and decomposed problems in OpenMDAO would be greatly appreciated.
In my code I used a top level assembly with a driver who's objective is to minimize the norm of c (the consistency constraints between targets and responses).

Assembly File

from openmdao.lib.datatypes.api import Float, Array
from openmdao.main.api import Assembly, set_as_top
from openmdao.lib.drivers.api import SLSQPdriver
import DisciplineTemplate
from numpy import ones, array, float

class Optimization0(Assembly):
    """my optimization template"""

    def configure(self):

        # Create Optimizer instance
        self.add('driver', SLSQPdriver())
        self.driver.accuracy = 1.0e-6

        self.add('temp', Float(1.0, iotype='in'))
        self.driver.add_parameter('temp', low=0., high=10.)

        # Create component instances
        self.add('Multipliers', DisciplineTemplate.Multipliers())

        # Iteration Hierarchy
        self.driver.workflow.add(['Multipliers'])

        # SLSQP Flags
        self.driver.iprint = 0

class Optimization1(Assembly):
    """my optimization template"""

    v = Array(array([1.0,1.0,1.0]), dtype=float, desc = "Penalty", iotype='in')
    w = Array(array([1.0,1.0,1.0]), dtype=float, desc = "Penalty", iotype='in')
    pR11 = Float(1.0, desc='Target', iotype='in')
    pR12 = Float(1.0, desc='Coupling', iotype='in')

    def configure(self):

        # Create Optimizer instance
        self.add('driver', SLSQPdriver())

        # Create component instances
        self.add('R01', DisciplineTemplate.R01())
        self.add('R02', DisciplineTemplate.R02())

        # Iteration Hierarchy
        self.driver.workflow.add(['R01', 'R02'])

        #x0[z4 z5 z7]
        self.driver.add_parameter('R01.z4', low=0., high=10.)
        self.driver.add_parameter(('R01.z5','R02.z5'), low = 0.0,   high=10.0)
        self.driver.add_parameter('R02.z7', low=0., high=10.)
        #T11=z3
        self.driver.add_parameter('R01.z3', low=0., high=10.)
        #T12=z6
        self.driver.add_parameter('R02.z6', low=0., high=10.)

        #Targets
        #self.create_passthrough('R01.z3', 'pT11')
        #self.create_passthrough('R02.z6', 'pT12')
        self.add('pT11', Float(1.0, iotype='in'))
        self.connect('R01.z3', 'pT11')
        self.add('pT12', Float(1.0, iotype='in'))
        self.connect('R02.z6', 'pT12')

        #Constraints
        self.driver.add_constraint('(R01.z3)**(-2) + (R01.z4)**(2) - (R01.z5)**(2) <= 0')
        self.driver.add_constraint('(R02.z5)**(2) + (R02.z6)**(-2) - (R02.z7)**(2) <= 0')

        # SLSQP Flags
        self.driver.iprint = 0

class Optimization2(Assembly):
    """my optimization template"""

    v = Array(array([1.0,1.0,1.0]), dtype=float, desc = "Penalty", iotype='in')
    w = Array(array([1.0,1.0,1.0]), dtype=float, desc = "Penalty", iotype='in')
    pT11 = Float(1.0, desc='Target', iotype='in')
    py12 = Float(1.0, desc='Coupling', iotype='in')

    def configure(self):

        # Create Optimizer instance
        self.add('driver', SLSQPdriver())

        # Create component instances
        self.add('R11', DisciplineTemplate.R11())

        # Iteration Hierarchy
        self.driver.workflow.add(['R11'])

        #x11[z8,z9,z10]
        self.driver.add_parameter('R11.z8', low=0., high=10.)
        self.driver.add_parameter('R11.z9', low=0., high=10.)
        self.driver.add_parameter('R11.z10', low=0., high=10.)
        #y11=z11
        self.driver.add_parameter('R11.z11', low=0., high=10.)

        #Out Targets
        #self.create_passthrough('R11.z3', 'pR11')
        self.add('pR11', Float(1.0, iotype='out'))
        self.connect('R11.z3', 'pR11')
        #self.create_passthrough('R11.z11', 'py11')
        self.add('py11', Float(1.0, iotype='out'))
        self.connect('R11.z11', 'py11')

        #Constraints
        self.driver.add_constraint('(R11.z8)**(2) + (R11.z9)**(2) - (R11.z11)**(2) <= 0')
        self.driver.add_constraint('(R11.z8)**(-2) + (R11.z10)**(2) - (R11.z11)**(2) <= 0')

        # SLSQP Flags
        self.driver.iprint = 0

class Optimization3(Assembly):
    """my optimization template"""

    v = Array(array([1.0,1.0,1.0]), dtype=float, desc = "Penalty", iotype='in')
    w = Array(array([1.0,1.0,1.0]), dtype=float, desc = "Penalty", iotype='in')
    pT12 = Float(1.0, desc='Target', iotype='in')
    py11 = Float(1.0, desc='Coupling', iotype='in')

    def configure(self):

        # Create Optimizer instance
        self.add('driver', SLSQPdriver())

        # Create component instances
        self.add('R12', DisciplineTemplate.R12())

        # Iteration Hierarchy
        self.driver.workflow.add(['R12'])

        #x12[z12,z13,z14]
        self.driver.add_parameter('R12.z12', low=0., high=10.)
        self.driver.add_parameter('R12.z13', low=0., high=10.)
        self.driver.add_parameter('R12.z14', low=0., high=10.)
        #y12=z11
        self.driver.add_parameter('R12.z11', low=0., high=10.)

        #Targets
        #self.create_passthrough('R12.z6', 'pR12')
        self.add('pR12', Float(1.0, iotype='out'))
        self.connect('R12.z6', 'pR12')
        #self.create_passthrough('R12.z11', 'py12')
        self.add('py12', Float(1.0, iotype='out'))
        self.connect('R12.z11', 'py12')

        #Constraints
        self.driver.add_constraint('(R12.z11)**(2) + (R12.z12)**(-2) - (R12.z13)**(2) <= 0')
        self.driver.add_constraint('(R12.z11)**(2) + (R12.z12)**(2) - (R12.z14)**(2) <= 0')

        # SLSQP Flags
        self.driver.iprint = 0

if __name__ == "__main__":

    #Hierarchy
    Loop = set_as_top(Optimization0())

    Loop.add('P0', Optimization1())
    Loop.driver.workflow.add('P0')

    Loop.add('P11', Optimization2())
    Loop.driver.workflow.add('P11')

    Loop.add('P12', Optimization3())
    Loop.driver.workflow.add('P12')

    #connections
    Loop.connect('[P11.pT11-P11.pR11,P12.pT12-P12.pR12,P12.py12-P11.py11]','Multipliers.c')

    Loop.connect('Multipliers.v','P0.v')
    Loop.connect('Multipliers.w','P0.w')
    Loop.connect('P11.pR11','P0.pR11')
    Loop.connect('P12.pR12','P0.pR12')

    Loop.connect('Multipliers.v','P11.v')
    Loop.connect('Multipliers.w','P11.w')
    Loop.connect('P0.pT11','P11.pT11')
    Loop.connect('P12.py12','P11.py12')

    Loop.connect('Multipliers.v','P12.v')
    Loop.connect('Multipliers.w','P12.w')
    Loop.connect('P0.pT12','P12.pT12')
    Loop.connect('P11.py11','P12.py11')

    #Objectives
    Loop.driver.add_objective('(P0.pR11-P0.pT11)**2 + (P0.pR12-P0.pT12)**2 + (P12.py12-P11.py11)**2')

    Loop.P0.driver.add_objective('(R01.z1-0.0)**2 + (R02.z2-0.0)**2 + v[0]*(R01.z3-pR11) + v[1]*(R02.z6-pR12) + w[0]*w[0]*(R01.z3-pR11)**2 + w[1]*w[1]*(R02.z6-pR12)**2')

    Loop.P11.driver.add_objective('v[0]*(pT11-R11.z3) + v[2]*(py12-R11.z11) + w[0]*w[0]*(pT11-R11.z3)**2 + w[2]*w[2]*(py12-R11.z11)**2')

    Loop.P12.driver.add_objective('v[1]*(pT12-R12.z6) + v[2]*(R12.z11-py11) + w[1]*w[1]*(pT12-R12.z6)**2 + w[2]*w[2]*(R12.z11-py11)**2')

    import time
    tt = time.time()

    Loop.run()

    official=[2.80,3.03,2.35,0.76,0.87,2.79,0.95,0.97,0.87,0.80,1.30,0.84,1.75,1.54]

    z=ones(14)
    z[0]=Loop.P0.R01.z1
    z[1]=Loop.P0.R02.z2
    z[2]=Loop.P0.R01.z3
    z[3]=Loop.P0.R01.z4
    z[4]=Loop.P0.R01.z5
    z[5]=Loop.P0.R02.z6
    z[6]=Loop.P0.R02.z7
    z[7]=Loop.P11.R11.z8
    z[8]=Loop.P11.R11.z9
    z[9]=Loop.P11.R11.z10
    z[10]=Loop.P11.R11.z11
    z[11]=Loop.P12.R12.z12
    z[12]=Loop.P12.R12.z13
    z[13]=Loop.P12.R12.z14

    for i in range(0, 14):
        print "\t\t z%d & %f & %f \\\\" % (i+1, z[i],official[i])

    print "\n"

    for i in range(0, 14):
        print "z%d=%f vs official value: %f" % (i+1, z[i], official[i])

    print "\nElapsed time: ", time.time()-tt, "seconds \n"

    print "Function calls R01: %d, R02: %d, R11: %d, R12: %d"%(Loop.P0.R01.exec_count, Loop.P0.R02.exec_count, Loop.P11.R11.exec_count, Loop.P12.R12.exec_count)

Components File

from openmdao.main.api import Component
from openmdao.lib.datatypes.api import Float, Array
from numpy import ones, array, float

class R01(Component):

    z3 = Float(1.0, iotype='in', desc='Coupling Design Variable')
    z4 = Float(1.0, iotype='in', desc='Local Design Variable')
    z5 = Float(1.0, iotype='in', desc='Shared Design Variable')

    z1 = Float(iotype='out', desc='Output of this Discipline')

    def execute(self):

        z3 = self.z3
        z4 = self.z4
        z5 = self.z5

        self.z1 = ( z3**2.0 + z4**(-2.0) + z5**2.0 )**(1.0/2.0)
        print "R01"
        print "z3=", z3
        print "z4=", z4
        print "z5=", z5        
        print "z1=", self.z1

class R02(Component):

    z5 = Float(1.0, iotype='in', desc='Shared Design Variable')
    z6 = Float(1.0, iotype='in', desc='Coupling Design Variable')
    z7 = Float(1.0, iotype='in', desc='Local Design Variable')

    z2 = Float(iotype='out', desc='Output of this Discipline')

    def execute(self):

        z5 = self.z5
        z6 = self.z6
        z7 = self.z7

        self.z2 = ( z5**2.0 + z6**(2.0) + z7**2.0 )**(1.0/2.0)
        print "R02"
        print "z5=", z5
        print "z6=", z6
        print "z7=", z7        
        print "z2=", self.z2

class R11(Component):

    z8 = Float(1.0, iotype='in', desc='Local Design Variable')
    z9 = Float(1.0, iotype='in', desc='Local Design Variable')
    z10 = Float(1.0, iotype='in', desc='Local Design Variable')
    z11 = Float(1.0, iotype='in', desc='Shared Design Variable')

    z3 = Float(iotype='out', desc='Output of this Discipline')

    def execute(self):

        z8 = self.z8
        z9 = self.z9
        z10 = self.z10
        z11 = self.z11

        self.z3 = ( z8**2.0 + z9**(-2.0) + z10**(-2.0) + z11**2.0 )**(1.0/2.0)
        print "R11"
        print "z8=", z8
        print "z9=", z9
        print "z10=", z10        
        print "z11=", z11
        print "z3=", self.z3

class R12(Component):

    z11 = Float(1.0, iotype='in', desc='Shared Design Variable')
    z12 = Float(1.0, iotype='in', desc='Local Design Variable')
    z13 = Float(1.0, iotype='in', desc='Local Design Variable')
    z14 = Float(1.0, iotype='in', desc='Local Design Variable')

    z6 = Float(iotype='out', desc='Output of this Discipline')

    def execute(self):

        z11 = self.z11
        z12 = self.z12
        z13 = self.z13
        z14 = self.z14

        self.z6 = ( z11**2.0 + z12**2.0 + z13**2.0 + z14**2.0 )**(1.0/2.0)
        print "R12"
        print "z11=", z11
        print "z12=", z12
        print "z13=", z13 
        print "z14=", z14        
        print "z6=", self.z6

class Multipliers(Component):
    c = Array(array([1.0,1.0,1.0]), dtype=float, desc = "Residuals", iotype='in')
    v = Array(array([1.0,1.0,1.0]), dtype=float, desc = "Penalty", iotype='out')
    w = Array(array([1.0,1.0,1.0]), dtype=float, desc = "Penalty", iotype='out')
    beta = Float(1.25, desc='beta')
    #c = [T11-R11 T12-R12 y12-y11]

    def execute(self):

        v = self.v
        w = self.w
        c = self.c
        beta = self.beta

        self.v = v + 2*w*w*c
        self.w = beta*w

        print "Multipliers"
        print "c=", self.c
        print "v=", self.v
        print "w=", self.w

asked 06 Aug '15, 17:24

Grasshopper's gravatar image

Grasshopper
3719


your ordering problem is related to your data passing. For instance you have

 Loop.connect('P0.pT11','P11.pT11')

and Loop.connect('P11.pR11','P0.pR11')

So you have a cycle between P0 and P1. OpenMDAO is probably ordering these in a different way then you would like. I think you want a jacobi like convergence behavior here, where each iteration just uses the values from the previous case and you loop around till it converges. You will need to look into using a the FixedPointIterator here. You can break one of your connections (the ones coming into P0) and then add the parameters and constraints to the FixedPointIterator accordingly. This should give you proper run order and actual convergence behavior.

You can look at our sellar_bliss example to see a good way to set this up.

link

answered 06 Aug '15, 20:00

justingray's gravatar image

justingray ♦♦
1.8k13

Follow this question

By Email:

Once you sign in you will be able to subscribe for any updates here

By RSS:

Answers

Answers and Comments

Markdown Basics

  • *italic* or _italic_
  • **bold** or __bold__
  • link:[text](http://url.com/ "title")
  • image?![alt text](/path/img.jpg "title")
  • numbered list: 1. Foo 2. Bar
  • to add a line break simply add two spaces to where you would like the new line to be.
  • basic HTML tags are also supported

Tags:

×39
×5
×3

Asked: 06 Aug '15, 17:24

Seen: 1,255 times

Last updated: 06 Aug '15, 20:00

powered by OSQA