add continuation on mma parameter for constraint relaxation

Ccore/onelab_optimize.py

@@ -130,7 +130,13 @@ def Optimize(Optimizer):
         conveks.mma.initialize(initialPoint, lowerBound, upperBound)
 
         # Set some options for MMA
+        if Optimizer == 1:
             conveks.mma.option.setNumber('General.Verbosity', 0) # def 4
+            conveks.mma.option.setNumber('SubProblem.isContinuation', 1)
+            conveks.mma.option.setNumber('SubProblem.cInitial'  , 1e2)
+            conveks.mma.option.setNumber('SubProblem.cFinal'    , 1e5)
+            conveks.mma.option.setNumber('SubProblem.cMultiply' , 10)
+            conveks.mma.option.setNumber('SubProblem.cIteration', 50)
         
         # Get iteration count (here it will be 1 - could be different in case of restart)
         it = conveks.mma.getOuterIteration()

Lbracket/topo.pro

@@ -9,6 +9,7 @@ DefineConstant [
   Optimizer = {0,Choices {0="Conveks - MMA",1="Conveks - GCMMA"},Name "Optimization/00Optimizer"}
   Opt_maxIter = {1600, Name "Optimization/01Max iterations"}
   Opt_maxChange = {0.01, Name "Optimization/02Max change"}
+  Opt_isInnerIter = {0, Name "Optimization/03Inner Iteration", Visible 0, Choices {0, 1}}
 
   // Density field parameters 
   densityFieldInit = {0.5, Name "Optimization/3Density/0Inital value"}
@@ -624,16 +625,16 @@ PostOperation {
         Format Table, StoreInVariable $vonMises_qp_pnorm, File "", Color "LightYellow",
         SendToServer StrCat[Opt_ResDir_Onelab,"pnorm of Von-Mises"] ] ;
 
-      If(Flag_PrintLevel > 2)
+      If(Flag_PrintLevel > 2 && !Opt_isInnerIter)
         CreateDir[Opt_ResDir];
         Print[ vonMises_qp, OnElementsOf Vol_Elast~{iP},
           File StrCat[Opt_ResDir,"vm_qp.pos"]];
       EndIf
-      If(Flag_PrintLevel > 4)
+      If(Flag_PrintLevel > 4 && !Opt_isInnerIter)
         Print[ vonMises, OnElementsOf Vol_Elast~{iP},
           File StrCat[Opt_ResDir,"vm.pos"]];
       EndIf
-      If(Flag_PrintLevel > 5)
+      If(Flag_PrintLevel > 5 && !Opt_isInnerIter)
         Print[ Young, OnElementsOf Vol_Elast~{iP},
           File StrCat[Opt_ResDir, "young.pos"] ];
         //Print[ Felast_pressure, OnElementsOf Sur_Force~{iP},
@@ -762,7 +763,7 @@ Resolution {
       InitSolution[LAM];
 
       // Show the density field
-      If(Flag_PrintLevel > 0)
+      If(Flag_PrintLevel > 0 && !Opt_isInnerIter)
         PostOperation[Get_DensityField];
       EndIf
 

Lbracket/topo.py

@@ -131,10 +131,18 @@ conveks.initialize()
 # Initialize the MMA optimizer
 conveks.mma.initialize(x, lowerBound, upperBound)
 
-
 # Set appropriate options for MMA
 conveks.mma.option.setNumber('General.Verbosity', 0)
-if Optimizer==0:conveks.mma.option.setNumber('SubProblem.move', 0.1)
+# if Optimizer==0:conveks.mma.option.setNumber('SubProblem.move', 0.1)
+# conveks.mma.option.setNumber('SubProblem.asymptotesRmax', 10)
+# conveks.mma.option.setNumber('SubProblem.asymptotesRmin', 0.01)
+# conveks.mma.option.setNumber('SubProblem.adaptSubproblem', 1)
+if Optimizer != 1:
+    conveks.mma.option.setNumber('SubProblem.isContinuation', 1)
+    conveks.mma.option.setNumber('SubProblem.cInitial'  , 1e2)
+    conveks.mma.option.setNumber('SubProblem.cFinal'    , 1e5)
+    conveks.mma.option.setNumber('SubProblem.cMultiply' , 10)
+    conveks.mma.option.setNumber('SubProblem.cIteration', 50)
 
 # Get iteration count (here it will be 1 - could be different in case of restart)
 it = conveks.mma.getOuterIteration()
@@ -154,17 +162,28 @@ c.setNumber('Optimization/Results/objective', value=objective)
 c.addNumberChoice('Optimization/Results/objective', value=objective) 
 c.setNumber('Optimization/Results/max(|Constraints|)', value=np.max(np.abs(constraints))) 
 c.addNumberChoice('Optimization/Results/max(|Constraints|)', value=np.max(np.abs(constraints))) 
+c.setNumber('Optimization/03Inner Iteration', value=0) 
+show_inner = 1
 
 # show current iteration
+# print 'iter.  inner iter.   obj.   max(constr.)   L2-norm(kkt)   point'
+# print '%3i  %3i   %2.4e   %2.4e   %2.4e'%(it-1, innerit, objective, np.max(np.abs(constraints)), kkt_norm),
+# print x[0:5]
 print 'iter.  inner iter.   obj.   max(constr.)   L2-norm(kkt)   point'
-print '%3i  %3i   %2.4e   %2.4e   %2.4e'%(it-1, innerit, objective, np.max(np.abs(constraints)), kkt_norm),
-print x[0:5]
+print '%3i  %3i  %.6e  %.6e  %.6e  %3i'%(it-1, innerit, objective, np.max(np.abs(constraints)), kkt_norm, 0),
+print '[',
+for xk in x[0:5]:
+    print '%.6e'%(float(xk)),
+print ']'
 
 while change > maxChange and it <= maxIter and c.getString('topo/Action') != 'stop':
 
     # Solve MMA
     conveks.mma.updateCurrentPoint(objective,grad_objective,constraints,grad_constraints)
 
+    # and get the total number of nonlinear loops
+    nlloops = conveks.mma.countNonlinearIterations()
+
     # Let us check now if the MMA approximation is conservative 
     # at the new point; and adapt the latter if it is not the case.
     if Optimizer == 1:
@@ -174,29 +193,42 @@ while change > maxChange and it <= maxIter and c.getString('topo/Action') != 'st
         conserv = conveks.mma.isConservativeSubProblem(obj_new, c_new) 
 
         innerit = 0
-        print '\t it.  new-point (xn)       f(xc)       f(xn)      max(c(xc))     max(c(xn))'
-        print '\t%3i'%(innerit),
-        print conveks.mma.getCurrentPoint()[0:5],
-        print '  %.4e  %.4e  %.4e  %.4e'%(objective,obj_new,np.max(constraints), np.max(c_new))
+        # print '\t it.  new-point (xn)       f(xc)       f(xn)      max(c(xc))     max(c(xn))'
+        # print '\t%3i'%(innerit),
+        # print conveks.mma.getCurrentPoint()[0:5],
+        # print '  %.4e  %.4e  %.4e  %.4e'%(objective,obj_new,np.max(constraints), np.max(c_new))
     
         while innerit <= 15 and conserv==0:
+            # Do not show fields (density, von-Mises, ...) during inner iterations
+            c.setNumber('Optimization/03Inner Iteration', value=1)
             innerit += 1
 
             # adapt the subproblem to make it conservative and generate a new point
             conveks.mma.updateCurrentPoint(objective,grad_objective,
                 objectiveAtUpdatedPoint=obj_new,constraintsAtUpdatedPoint=c_new)
 
+            # and get the total number of nonlinear loops
+            nlloops = conveks.mma.countNonlinearIterations()
+
             # evaluate the objective and the constraints at the new point
             # and check if the subproblem is conservative at this point
             obj_new,c_new = get_objective(conveks.mma.getCurrentPoint())
             conserv = conveks.mma.isConservativeSubProblem(obj_new, c_new) 
 
             # show inner iteration
-            print '\t%3i'%(innerit),
-            print conveks.mma.getCurrentPoint()[0:5],
-            print '  %.4e  %.4e  %.4e  %.4e'%(objective,obj_new,np.max(constraints), np.max(c_new))
+            # print '\t%3i'%(innerit),
+            # print conveks.mma.getCurrentPoint()[0:5],
+            # print '  %.4e  %.4e  %.4e  %.4e  %3i'%(objective,obj_new,np.max(constraints), np.max(c_new), nlloops)
+            if show_inner:
+                print '%3i  %3i '%(it, innerit),
+                print '%.6e  %.6e  %.6e  %3i'%(obj_new,np.max(np.abs(c_new)), kkt_norm, nlloops),
+                print '[',
+                for xk in conveks.mma.getCurrentPoint()[0:5]:
+                    print '%.6e'%(float(xk)),
+                print ']'
 
         inneritTot += innerit
+        c.setNumber('Optimization/03Inner Iteration', value=0) 
 
     # get (copy of) current point
     x = conveks.mma.getCurrentPoint()
@@ -210,8 +242,13 @@ while change > maxChange and it <= maxIter and c.getString('topo/Action') != 'st
     kkt_norm = conveks.mma.getKKTNorm(constraints,grad_objective,grad_constraints)
     
     # show the current state
-    print '%3i  %3i  %2.4e  %2.4e  %2.4e'%(it, innerit, objective, np.max(np.abs(constraints)), kkt_norm),
-    print x[0:5]
+    # print '%3i  %3i  %2.4e  %2.4e  %2.4e  %4i'%(it, innerit, objective, np.max(np.abs(constraints)), kkt_norm, nlloops),
+    # print x[0:5]
+    print '%3i  %3i  %.6e  %.6e  %.6e  %3i'%(it, innerit, objective, np.max(np.abs(constraints)), kkt_norm, nlloops),
+    print '[',
+    for xk in x[0:5]:
+        print '%.6e'%(float(xk)),
+    print ']'
 
     # get the outer iteration count here as well as the 
     # change between two successive points 

Team25/shape.py

@@ -155,6 +155,13 @@ conveks.mma.initialize(initialPoint, lowerBound, upperBound)
 
 # Set some options for MMA
 conveks.mma.option.setNumber('General.Verbosity', 0*4)
+if Optimizer != 1:
+    conveks.mma.option.setNumber('General.Verbosity', 0) # def 4
+    conveks.mma.option.setNumber('SubProblem.isContinuation', 1)
+    conveks.mma.option.setNumber('SubProblem.cInitial'  , 1e2)
+    conveks.mma.option.setNumber('SubProblem.cFinal'    , 1e5)
+    conveks.mma.option.setNumber('SubProblem.cMultiply' , 10)
+    conveks.mma.option.setNumber('SubProblem.cIteration', 50)
 
 # Get iteration count (here it will be 1 - could be different in case of restart)
 it = conveks.mma.getOuterIteration()