From d431fb3c6cd520569e4638daf73100173b21a1b1 Mon Sep 17 00:00:00 2001
From: ErinKuci <Erin.Kuci@ulg.ac.be>
Date: Sat, 12 Jan 2019 15:20:01 +0100
Subject: [PATCH] add the conservative approximation of MMA
---
Ccore/ccore_common.pro | 3 +-
Ccore/onelab_optimize.py | 48 +++++++-
Lbracket/topo.pro | 40 +++++--
Lbracket/topo.py | 197 ++++++++++++++++++++++----------
Team25/shape.geo | 10 +-
Team25/shape.pro | 9 ++
Team25/shape.py | 240 ++++++++++++++++++++++++++++-----------
7 files changed, 396 insertions(+), 151 deletions(-)
diff --git a/Ccore/ccore_common.pro b/Ccore/ccore_common.pro
index 6a69180..c6d1186 100644
--- a/Ccore/ccore_common.pro
+++ b/Ccore/ccore_common.pro
@@ -61,7 +61,8 @@ DesignSpace =
Optimizer =
DefineNumber[0, Name "Optimization/Algorithm",
Choices {0="scipy",
- 1="conveks"}];
+ 1="conveks-MMA",
+ 2="conveks-GCMMA"}];
w =
DefineNumber[0, Name "Optimization/Results/w",
Graph "030000000000000000000000000000000000"];
diff --git a/Ccore/onelab_optimize.py b/Ccore/onelab_optimize.py
index e15cc02..6833da8 100644
--- a/Ccore/onelab_optimize.py
+++ b/Ccore/onelab_optimize.py
@@ -112,7 +112,7 @@ def Optimize(Optimizer):
print 'Stopped by Gmsh'
return
print res
- elif Optimizer==1: # 'mma'
+ elif Optimizer==1 or Optimizer==2: # 'mma'
numVariables = len(x)
initialPoint = np.zeros(numVariables)
lowerBound = np.zeros(numVariables)
@@ -134,7 +134,10 @@ def Optimize(Optimizer):
# Get iteration count (here it will be 1 - could be different in case of restart)
it = conveks.mma.getOuterIteration()
- change_mma = 1.0
+ innerit = 0; inneritTot = 0; change_mma = 1.0
+ all_objectives = []
+ print 'outer iter. inner iter. change obj. constr.'
+
while it <= maxIter :
c.setNumber('Optimization/01Current iteration', value=it, readOnly=1,
attributes={'Highlight':'LightYellow'})
@@ -145,7 +148,12 @@ def Optimize(Optimizer):
obj = fobj(xFromMMA)
if (c.getString(c.name+'/Action') == 'stop'): break
grad_obj = grad(xFromMMA)
-
+ constraints = np.array([np.sum(xFromMMA)/np.sum(upperBound)-1.0])
+ grad_constraints = np.ones(numVariables)/np.sum(upperBound)
+ all_objectives.append(obj)
+ #print '%3i %3i %.4e %.4e %.4e'%(it,innerit,change_mma,obj,constraints[0])
+ #print xFromMMA
+
if it == 1:
fscale = 1.0
# fscale = 1.0 / np.abs(obj)
@@ -164,7 +172,37 @@ def Optimize(Optimizer):
# call MMA and update the current point
#conveks.mma.setMoveLimits(lowerBound, upperBound, 0.01) # parametre interessant
- conveks.mma.updateCurrentPoint(grad_objective)
+ #conveks.mma.updateCurrentPoint(grad_objective)
+ conveks.mma.updateCurrentPoint(objective,grad_objective,constraints,grad_constraints)
+
+ if Optimizer==2: # we use GCMMA, so update the approximation to be conservative
+ # Evaluate the objective at the new point
+ xFromMMA = conveks.mma.getCurrentPoint()
+ objective_new = fobj(xFromMMA)
+ constraints_new = np.array([np.sum(xFromMMA)/np.sum(upperBound)-1.0])
+
+ # check if the approximation is conservative at the current point
+ conserv = conveks.mma.isConservativeSubProblem(objective_new,constraints_new)
+
+ innerit = 0
+ print '\t',innerit,xFromMMA,objective_new,objective
+ while innerit <= 15 and conserv==0:
+ innerit += 1
+
+ # adapt the subproblem to make it conservative and generate a new point
+ conveks.mma.updateCurrentPoint(objective,grad_objective,constraints,grad_constraints,
+ objectiveAtUpdatedPoint=objective_new,constraintsAtUpdatedPoint=constraints_new)
+
+ # evaluate the objective and the constraints at the new point
+ xFromMMA = conveks.mma.getCurrentPoint()
+ objective_new = fobj(xFromMMA)
+ constraints_new = np.array([np.sum(xFromMMA)/np.sum(upperBound)-1.0])
+ print '\t',innerit,xFromMMA,objective_new,objective,constraints_new,constraints
+
+ # and check if the subproblem is conservative at this point
+ conserv = conveks.mma.isConservativeSubProblem(objective_new,constraints_new)
+ inneritTot += innerit
+
change_mma = conveks.mma.getDesignChange()
if change_mma < maxChange: break
it = conveks.mma.getOuterIteration()
@@ -172,6 +210,8 @@ def Optimize(Optimizer):
print
print 'Converged: %4d %3.6f %3.6f %3.6f' %(it, obj, change_rel, change_mma),
print xFromMMA, grad_obj
+ print 'objective'
+ for i in all_objectives:print i
# This should be called at the end of MMA
conveks.mma.finalize()
diff --git a/Lbracket/topo.pro b/Lbracket/topo.pro
index c16cfe3..9aa2aea 100644
--- a/Lbracket/topo.pro
+++ b/Lbracket/topo.pro
@@ -3,9 +3,19 @@ Struct OPT::YOUNG_LAW [ Enum, simp, modifiedSimp, ramp, none];
DefineConstant [
Opt_ResDir = StrCat[CurrentDir,"res_opt/"]
Opt_ResDir_Onelab = "Optimization/Results/"
- densityFieldInit = {0.5, Name "Optimization/3Density/0Inital value"}
- Opt_filter_radius = {0.008, Name "Optimization/3Density/2filter radius"}
Flag_PrintLevel = {1, Name "General/Verbosity", Visible 1}
+
+ // Optimization parameters
+ 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"}
+
+ // Density field parameters
+ densityFieldInit = {0.5, Name "Optimization/3Density/0Inital value"}
+ Opt_filter_exists = {1, Name "Optimization/3Density/1filter", Choices{0,1}}
+ Opt_filter_radius = {0.008, Name "Optimization/3Density/2filter radius", Visible Opt_filter_exists}
+
+ // Interpolation parameters
Flag_opt_matlaw = {OPT::YOUNG_LAW.modifiedSimp,
Choices {
OPT::YOUNG_LAW.simp="simp",
@@ -19,12 +29,21 @@ DefineConstant [
|| Flag_opt_matlaw == OPT::YOUNG_LAW.ramp
|| Flag_opt_matlaw == OPT::YOUNG_LAW.modifiedSimp) }
Flag_opt_stress_penal = {2.5,
- Name "Optimization/5Young Interpolation/2penalization of stress",
+ Name "Optimization/6Von-Mises parameters/1penalization of stress",
Visible Flag_opt_matlaw != OPT::YOUNG_LAW.none}
+ PARAM_pnorm = {8,Name "Optimization/6Von-Mises parameters/2Pnorm penalization"}
+
+ // Model parameters
Flag_EPC = 1
iP = 1
+ Flag_getGradient = 1
];
+fobj = DefineNumber[0, Name "Optimization/Results/objective",
+ Graph "0300000000000000000000000000000000000"];
+maxConstr = DefineNumber[0, Name "Optimization/Results/max(|Constraints|)",
+ Graph "0000000003000000000000000000000000000"];
+
Group {
Sur_Clamp~{iP} = Region[{1001}];
Vol_Force~{iP} = Region[{}];
@@ -535,11 +554,6 @@ Return
// We also make use of the adjoint method to compute their
// sensitivity wit respect to the density field.
-DefineConstant[
- PARAM_pnorm = {8,
- Name "Optimization/7Parameter(s) of performances/Von-Mises/Pnorm penalization"}
-];
-
PostProcessing{
{ Name ObjectiveConstraints; NameOfFormulation Elast_u;
Quantity {
@@ -763,10 +777,12 @@ Resolution {
PostOperation[Get_ObjectiveConstraints];
// Compute the sensitivity of objective/constraints
- PostOperation[Get_GradOf_Volume];
- Generate[LAM];
- SolveAgainWithOther[LAM,SYS];
- PostOperation[Get_GradOf_VonMises];
+ If(Flag_getGradient)
+ PostOperation[Get_GradOf_Volume];
+ Generate[LAM];
+ SolveAgainWithOther[LAM,SYS];
+ PostOperation[Get_GradOf_VonMises];
+ EndIf
}
}
}
diff --git a/Lbracket/topo.py b/Lbracket/topo.py
index 1606aae..bc333d1 100644
--- a/Lbracket/topo.py
+++ b/Lbracket/topo.py
@@ -39,15 +39,15 @@ c.openProject(topo_geo)
# dry getdp run (without -solve or -pos option) to get model parameters in the GUI
c.runSubClient('myGetDP', mygetdp)
-# define now optimization parameters
+# get now optimization parameters
# some of them as Onelab parameter, to be editable in the GUI
-maxIter = c.defineNumber('Optimization/00Max iterations', value=1600)
-maxChange = c.defineNumber('Optimization/01Max change', value=1e-2)
-filtering = c.defineNumber('Optimization/02Filtering/00Filter densities', value=1, choices=[0,1])
+Optimizer = c.getNumber('Optimization/00Optimizer')
+maxIter = c.getNumber('Optimization/01Max iterations')
+maxChange = c.getNumber('Optimization/02Max change')
+filtering = c.getNumber('Optimization/3Density/1filter')
# end of check phase. We're done if we do not start the optimization
-if c.action == 'check' :
- exit(0)
+if c.action == 'check': exit(0)
# mesh the geometry to create the density (background) mesh
c.runSubClient('myGmsh', mygmsh + '-2')
@@ -78,6 +78,44 @@ def setElementTable(var, ele, val):
data = np.insert(data, 0, numElements)
c.setNumber(var, values=data, visible=0)
+def get_objective(xFromOpti, flag_grad=1):
+ if filtering:
+ setElementTable('Optimization/Results/filterInput', elements, xFromOpti)
+ c.runSubClient('myGetDP', mygetdp + '-solve FilterDensity')
+ xFromOpti = np.array(c.getNumbers('Optimization/Results/filterOutput')[2::2])
+
+ setElementTable('Optimization/Results/designVariable', elements, xFromOpti)
+
+ c.runSubClient('myGetDP', mygetdp + '-solve ComputePerformancesAndGradient '\
+ + '-setnumber Flag_getGradient '+str(flag_grad))
+
+ # get the value of the objective function and of the constraints
+ # as well as their sensitivity with respect to design variables at `x'
+ objective = c.getNumber('Optimization/Results/pnorm of Von-Mises')
+ constraints = np.array([c.getNumber('Optimization/Results/VolumeDomain')\
+ / (0.5*c.getNumber('Optimization/Results/VolumeDomainFull')) - 1.0])
+
+ return objective, constraints
+
+def get_grad(xFromOpti):
+ grad_objective = np.asarray(c.getNumbers('Optimization/Results/GradOfVonMises')[2::2])
+ grad_constraints = np.asarray(c.getNumbers('Optimization/Results/GradOfVolume')[2::2])\
+ /(0.5*c.getNumber('Optimization/Results/VolumeDomainFull'))
+ if filtering:
+ # apply the chain rule for the sensitivity of objective
+ setElementTable('Optimization/Results/filterInput', elements, grad_objective)
+ c.runSubClient('myGetDP', mygetdp + '-solve FilterDensity')
+ grad_objective = np.array(c.getNumbers('Optimization/Results/filterOutput')[2::2])
+
+ # apply the chain rule for the sensitivity of constraints
+ d_con_dx = grad_constraints.reshape((1, numVariables))
+ for k,dfe in enumerate(d_con_dx):
+ setElementTable('Optimization/Results/filterInput', elements, dfe)
+ c.runSubClient('myGetDP', mygetdp + '-solve FilterDensity')
+ d_con_dx[k] = np.array(c.getNumbers('Optimization/Results/filterOutput')[2::2])
+ grad_constraints = d_con_dx.reshape((1 * numVariables,))
+ return grad_objective, grad_constraints
+
# number of design variables and number of constraints
numVariables = len(x)
@@ -93,77 +131,116 @@ conveks.initialize()
# Initialize the MMA optimizer
conveks.mma.initialize(x, lowerBound, upperBound)
+
# Set appropriate options for MMA
-conveks.mma.option.setNumber('General.Verbosity', 4)
-conveks.mma.option.setNumber('SubProblem.move', 0.1)
+conveks.mma.option.setNumber('General.Verbosity', 0)
+if Optimizer==0:conveks.mma.option.setNumber('SubProblem.move', 0.1)
# Get iteration count (here it will be 1 - could be different in case of restart)
it = conveks.mma.getOuterIteration()
-change = 1.0
+innerit = 0; inneritTot = 0; change = 1.0; kkt_norm=1.0
+all_objectives = []
-while change > maxChange and it <= maxIter and c.getString('topo/Action') != 'stop':
+# get (copy of) current point
+x = conveks.mma.getCurrentPoint()
- c.setNumber('Optimization/01Current iteration', value=it, readOnly=1,
- attributes={'Highlight':'LightYellow'})
- c.setNumber('Optimization/02Current change', value=change, readOnly=1,
- attributes={'Highlight':'LightYellow'})
+# compute objective function and constraints
+# as well as their sensitivity with respect to design variables at `x'
+objective,constraints = get_objective(x)
+grad_objective,grad_constraints = get_grad(x)
- # get (copy of) current point
- x = conveks.mma.getCurrentPoint()
+# send data to onelab server
+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)))
- if filtering:
- setElementTable('Optimization/Results/filterInput', elements, x)
- c.runSubClient('myGetDP', mygetdp + '-solve FilterDensity')
- x = np.array(c.getNumbers('Optimization/Results/filterOutput')[2::2])
+# 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]
- setElementTable('Optimization/Results/designVariable', elements, x)
+while change > maxChange and it <= maxIter and c.getString('topo/Action') != 'stop':
- # compute objective function and constraints as well as their sensitivity
- # with respect to design variables at `x'
- c.runSubClient('myGetDP', mygetdp + '-solve ComputePerformancesAndGradient')
+ # Solve MMA
+ conveks.mma.updateCurrentPoint(objective,grad_objective,constraints,grad_constraints)
+
+ # 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:
+ # Evaluate the objective at the new point and check if
+ # the approximation is conservative at the current point
+ obj_new,c_new = get_objective(conveks.mma.getCurrentPoint())
+ 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))
+
+ while innerit <= 15 and conserv==0:
+ 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)
+
+ # 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))
+
+ inneritTot += innerit
- # get the value of the objective function and of the constraints
+ # get (copy of) current point
+ x = conveks.mma.getCurrentPoint()
+
+ # compute objective function and constraints
# as well as their sensitivity with respect to design variables at `x'
- objective = c.getNumber('Optimization/Results/pnorm of Von-Mises')
- constraints = np.array([c.getNumber('Optimization/Results/VolumeDomain')\
- / (0.5*c.getNumber('Optimization/Results/VolumeDomainFull')) - 1.0])
- grad_objective = np.array(c.getNumbers('Optimization/Results/GradOfVonMises')[2::2])
- grad_constraints = np.asarray(c.getNumbers('Optimization/Results/GradOfVolume')[2::2])\
- /(0.5*c.getNumber('Optimization/Results/VolumeDomainFull'))
+ objective,constraints = get_objective(x)
+ grad_objective,grad_constraints = get_grad(x)
- if filtering:
- # apply the chain rule for the sensitivity of objective
- setElementTable('Optimization/Results/filterInput', elements, grad_objective)
- c.runSubClient('myGetDP', mygetdp + '-solve FilterDensity')
- grad_objective = np.array(c.getNumbers('Optimization/Results/filterOutput')[2::2])
+ # evaluate the L2-norm of KKT
+ kkt_norm = conveks.mma.getKKTNorm(constraints,grad_objective,grad_constraints)
- # apply the chain rule for the sensitivity of constraints
- d_con_dx = grad_constraints.reshape((len(constraints), numVariables))
- for k,dfe in enumerate(d_con_dx):
- setElementTable('Optimization/Results/filterInput', elements, dfe)
- c.runSubClient('myGetDP', mygetdp + '-solve FilterDensity')
- d_con_dx[k] = np.array(c.getNumbers('Optimization/Results/filterOutput')[2::2])
- grad_constraints = d_con_dx.reshape((len(constraints) * numVariables,))
-
- if it == 1: fscale = 1.0 / np.abs(objective)
- objective *= fscale
- grad_objective *= fscale
-
- print('*'*50)
- print('iteration: ', it)
- print('change: ', change)
- print('objective:', objective)
- print('constraints:', constraints)
- c.sendInfo('Optimization: it. {}, obj. {}, constr. {}, change {}'.format(it,objective,constraints[0], change))
-
- # call MMA and update the current point
- #conveks.mma.setMoveLimits(lowerBound, upperBound, 0.1)
- conveks.mma.updateCurrentPoint(grad_objective, constraints, grad_constraints)
- change = conveks.mma.getDesignChange()
+ # 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]
+
+ # get the outer iteration count here as well as the
+ # change between two successive points
it = conveks.mma.getOuterIteration()
+ change = conveks.mma.getDesignChange()
+
+ # send data to the onelab server
+ c.setNumber('Optimization/Results/01Current iteration', value=it, readOnly=1)
+ c.setNumber('Optimization/Results/02Current change', value=change, readOnly=1)
+ c.setNumber('Optimization/Results/03KKT L2-norm', value=kkt_norm, readOnly=1)
+ 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)))
+
+# show the problem at the converged state
+print '*'*80
+print 'Optimization has converged after %i iterations.'%(it-1)
+print ' > point :',x[0:5]
+print ' > objective : %.4e'%(objective)
+print ' > max(constraints): %.4e'%(np.max(np.abs(constraints)))
+print ' > L2-norm of kkt : %.4e'%(kkt_norm)
+print ' > inner iterations: %i'%(inneritTot)
+print '*'*80
# This should be called at the end of MMA
conveks.mma.finalize()
# This should be called at the end
-conveks.finalize()
\ No newline at end of file
+conveks.finalize()
+
+print all_objectives
\ No newline at end of file
diff --git a/Team25/shape.geo b/Team25/shape.geo
index f20bc36..097ccdc 100644
--- a/Team25/shape.geo
+++ b/Team25/shape.geo
@@ -4,7 +4,7 @@ DefineConstant[
R1 = 6*mm
L2 = 18*mm
L3 = 15*mm
- L4 = {L2/L3*Sqrt[L3^2-(12.5*mm)^2], Name "Constructive parameters/L4"}
+ L4 = {L2/L3*Sqrt[L3^2-(12.5*mm)^2], Name "Optimization/parameters/L4", Closed 1}
angleMeasure = 50*deg
lc = 12*mm
lcd = lc/15
@@ -20,7 +20,7 @@ Point(1) = {0, 0, 0, lcd};
For k In {0:10}
theta = k * 90/10 * deg;
- DefineConstant[Rs~{k} = {R1, Name Sprintf["Constructive parameters/spline radius %g",k+1]}];
+ DefineConstant[Rs~{k} = {R1, Name Sprintf["Optimization/parameters/spline radius %g",k+1], Closed 1}];
Point(2020+k) = {Rs~{k}*Cos[theta], Rs~{k}*Sin[theta], 0, lcd};
EndFor
@@ -33,7 +33,7 @@ BSpline(4022) = {2026, 2027, 2028, 2029, 2030};
For k In {0:9}
theta = k * 50/10 * deg;
DefineConstant[Rso~{k} = {L2*L3/Sqrt[(L3*Cos[theta])^2+(L2*Sin[theta])^2],
- Name Sprintf["Constructive parameters/outer spline radius %g",k+1]}];
+ Name Sprintf["Optimization/parameters/outer spline radius %g",k+1], Closed 1}];
Point(6020+k) = {Rso~{k}*Cos[theta], Rso~{k}*Sin[theta], 0, lcd};
EndFor
@@ -64,8 +64,8 @@ Point(1022) = {(9.5+2.25)*mm, 0, 0, lcd};
Point(23) = {20*mm, 12.5*mm, 0, lcd};
-//DefineConstant[edge_x_1 = {20*mm-L4, Name "Constructive parameters/edge radius 1"}];
-//DefineConstant[edge_y_1 = {(12.5-2)*mm, Name "Constructive parameters/edge radius 2"}];
+//DefineConstant[edge_x_1 = {20*mm-L4, Name "Optimization/parameters/edge radius 1"}];
+//DefineConstant[edge_y_1 = {(12.5-2)*mm, Name "Optimization/parameters/edge radius 2"}];
//Point(25) = {edge_x_1, edge_y_1, 0, lcd};
Circle(1016) = {1022, 1, 1021};
diff --git a/Team25/shape.pro b/Team25/shape.pro
index d27346f..8cdb302 100644
--- a/Team25/shape.pro
+++ b/Team25/shape.pro
@@ -73,8 +73,17 @@ DefineConstant [
Nb_max_iter = 30
relaxation_factor = 1
stop_criterion = 1e-5
+ Optimizer = {0,Choices {0="Conveks - MMA",1="Conveks - GCMMA"},Name "Optimization/00Optimizer"}
+ Opt_maxIter = {1600, Name "Optimization/01Max iterations"}
+ Opt_maxChange = {0.0001, Name "Optimization/02Max change"}
+ Opt_maxKKT = {0.001, Name "Optimization/03Max KKT norm"}
];
+fobj = DefineNumber[0, Name "Optimization/Results/objective",
+ Graph "0300000000000000000000000000000000000"];
+maxConstr = DefineNumber[0, Name "Optimization/Results/max(|Constraints|)",
+ Graph "0000000003000000000000000000000000000"];
+
Group {
// One starts by giving explicit meaningful names to the Physical regions
// defined in the "shape.geo" mesh file.
diff --git a/Team25/shape.py b/Team25/shape.py
index 310805d..60b3200 100644
--- a/Team25/shape.py
+++ b/Team25/shape.py
@@ -21,13 +21,13 @@ if(not len(mygetdp)):
c.sendError('Please run a GetDP model interactively once with Gmsh to ' +
'initialize the solver location')
exit(0)
-
c.sendInfo('Will use gmsh={0} and getdp={1}'.format(mygmsh, mygetdp))
# append correct path to model file names
-file_geo = c.getPath('shape.geo')
-file_msh = c.getPath('shape.msh')
-file_pro = c.getPath('shape.pro')
+modelName = 'shape'
+file_geo = c.getPath(modelName+'.geo')
+file_msh = c.getPath(modelName+'.msh')
+file_pro = c.getPath(modelName+'.pro')
# build command strings with appropriate parameters and options
mygetdp += ' -p 0 ' + file_pro + ' -msh ' + file_msh + ' '
@@ -39,56 +39,50 @@ c.openProject(file_geo)
# dry getdp run (without -solve or -pos option) to get model parameters in the GUI
c.runSubClient('myGetDP', mygetdp)
-# define optimization parameters as Onelab parameter (editable in the GUI)
-maxIter = c.defineNumber('Optimization/00Max iterations', value=100)
-maxChange = c.defineNumber('Optimization/01Max change', value=1e-5)
+# read optimization parameters as Onelab parameter (editable in the GUI)
+Optimizer = c.getNumber('Optimization/00Optimizer')
+maxIter = c.getNumber('Optimization/01Max iterations')
+maxChange = c.getNumber('Optimization/02Max change')
+maxKKTNorm = c.getNumber('Optimization/03Max KKT norm')
# end of check phase. We're done if we do not start the optimization
-if c.action == 'check' :
- exit(0)
+if c.action == 'check' :exit(0)
-x = {}
+# Let us define now the design variables as a list
+# for which the entries are defiend as follows:
+# ['onelab/path', value, lower-bound, upper-bound]
+allParameters = []
+xe = 6.0e-3
for k in xrange(11):
- xe = 6.0e-3
- x[k] = ['Constructive parameters/spline radius '+str(k+1),xe, xe-4.0e-3, xe+4.0e-3, 'Rs_'+str(k)]
-sizeX = len(x.keys())
+ allParameters.append(['Optimization/parameters/spline radius '+str(k+1),xe, xe-4.0e-3, xe+4.0e-3, 'Rs_'+str(k)])
L2 = 18.0e-3;L3 = 15.0e-3
for k in xrange(10):
theta = float(k) * 50.0/10.0 * np.pi/180.0
xe = L2*L3/((L3*np.cos(theta))**2.0+(L2*np.sin(theta))**2.0)**0.5
- x[sizeX+k] = ['Constructive parameters/outer spline radius '+str(k+1), xe, xe-4.0e-3, xe+1.5e-3, 'Rso_'+str(k)]
+ allParameters.append(['Optimization/parameters/outer spline radius '+str(k+1), xe, xe-4.0e-3, xe+1.5e-3, 'Rso_'+str(k)])
-x[len(x.keys())] = ['Constructive parameters/L4', 0.00994987, 0.00994987-4.0e-3, xe+4.0e-3, 'L4']
+allParameters.append(['Optimization/parameters/L4', 0.00994987, 0.00994987-4.0e-3, xe+4.0e-3, 'L4'])
+
+# we create a dictionnary based on the list of design variables
+x = {}
+for k in xrange(len(allParameters)):
+ x[k] = allParameters[k]
+# some utility functions:
+# (1) to read the values from a table
def readSimpleTable(path):
with open(path) as FileObj:
return np.array([float(lines.split()[-1]) for lines in FileObj])
+# (2) to set the values of a node table
def setNodeTable(var, nodes, val):
data = np.ravel(np.column_stack((nodes, val)))
data = np.insert(data, 0, float(len(nodes)))
c.setNumber(var, values=data, visible=0)
-def getVelocityField(xvar, perturb=1e-6):
- for id, var in xvar.iteritems():
- c.runNonBlockingSubClient('myGmsh', mygmsh \
- + ' -setnumber PerturbMesh 1'\
- + ' -setnumber VelocityTag '+str(id)\
- + ' -setnumber '+var[-1]+' '+str(var[1]+perturb)\
- + ' -run')
- c.waitOnSubClients()
- # send the x,y,z components of each velocity field
- for label, var in x.iteritems():
- d = c.getNumbers('Optimization/Results/velocity_'+str(label))
- if label==0:nodes = np.array(d[1::4])
- setNodeTable('Optimization/Results/velocity_'+str(label)+'_1',nodes,np.array(d[2::4]))
- setNodeTable('Optimization/Results/velocity_'+str(label)+'_2',nodes,np.array(d[3::4]))
- setNodeTable('Optimization/Results/velocity_'+str(label)+'_3',nodes,np.array(d[4::4]))
- c.setNumber(var[0],value=x[label][1])
-
# number of design variables and number of constraints
-numVariables = len(x.keys())
+numVariables = len(x)
# initial value of the variables in the design space,
# the lower bound for design variables,
@@ -97,34 +91,16 @@ initialPoint = np.zeros(numVariables)
lowerBound = np.zeros(numVariables)
upperBound = np.zeros(numVariables)
for label, var in x.iteritems():
- initialPoint[label] = var[1]
- lowerBound[label] = var[2]
- upperBound[label] = var[3]
-
-# Initialize conveks
-conveks.initialize()
-
-# Initialize the MMA optimizer
-conveks.mma.initialize(initialPoint, lowerBound, upperBound)
-
-# Set some options for MMA
-conveks.mma.option.setNumber('General.Verbosity', 4)
-
-# Get iteration count (here it will be 1 - could be different in case of restart)
-it = conveks.mma.getOuterIteration()
-change = 1.0
-
-while it <= maxIter and c.getString('shape/Action') != 'stop':
-
- c.setNumber('Optimization/01Current iteration', value=it, readOnly=1,
- attributes={'Highlight':'LightYellow'})
-
- # get (copy of) current point
- xFromMMA = conveks.mma.getCurrentPoint()
+ initialPoint[label] = var[1]
+ lowerBound[label] = var[2]
+ upperBound[label] = var[3]
+# We define the function which returns the objective and
+# constraints values (as a list) for a given point 'xFromOpti'.
+def get_objective(xFromOpti):
# send the current point to GetDP model
for label, var in x.iteritems():
- x[label][1] = xFromMMA[label]
+ x[label][1] = xFromOpti[label]
c.setNumber(var[0],value=x[label][1])
# reload the geometry in GUI to see the geometrical changes
@@ -133,12 +109,31 @@ while it <= maxIter and c.getString('shape/Action') != 'stop':
# mesh the geometry
c.runSubClient('myGmsh', mygmsh + ' -2')
- # compute objective function and constraints
+ # run Fem and evaluate the objective function
c.runSubClient('myGetDP', mygetdp + '-solve GetPerformances')
objective = np.sum(readSimpleTable(c.getPath('res_opt/w.txt')))
+ constraints = np.array([np.sum(xFromOpti)/np.sum(upperBound)-1]);
+ return objective,constraints
+# We also define the function which returns the gradient of the
+# objective and constraints (as a C-style array) for a given point.
+def get_grad(perturb=1e-6):
# generate the velocity field of each design variable
- getVelocityField(x)
+ for id, var in x.iteritems():
+ c.runNonBlockingSubClient('myGmsh', mygmsh \
+ + ' -setnumber PerturbMesh 1'\
+ + ' -setnumber VelocityTag '+str(id)\
+ + ' -setnumber '+var[-1]+' '+str(var[1]+perturb)\
+ + ' -run')
+ c.waitOnSubClients()
+ # send the x,y,z components of each velocity field
+ for label, var in x.iteritems():
+ d = c.getNumbers('Optimization/Results/velocity_'+str(label))
+ if label==0:nodes = np.array(d[1::4])
+ setNodeTable('Optimization/Results/velocity_'+str(label)+'_1',nodes,np.array(d[2::4]))
+ setNodeTable('Optimization/Results/velocity_'+str(label)+'_2',nodes,np.array(d[3::4]))
+ setNodeTable('Optimization/Results/velocity_'+str(label)+'_3',nodes,np.array(d[4::4]))
+ c.setNumber(var[0],value=x[label][1])
# compute the sensitivity with respect to design variables at `x'
for dv, var in x.iteritems():
@@ -148,20 +143,127 @@ while it <= maxIter and c.getString('shape/Action') != 'stop':
c.waitOnSubClients()
grad_objective = np.asarray([np.sum(readSimpleTable(c.getPath('res_opt/Grad_w_wrt_dv_'+str(dv)+'.txt')))\
for dv in xrange(numVariables)])
+ grad_constraints = np.ones(numVariables)/np.sum(upperBound);
+
+ return grad_objective,grad_constraints
- print('*'*50)
- print('iteration: ', it)
- print('change: ', change)
- print('current point:', xFromMMA)
- print('objective:', objective)
- c.sendInfo('Optimization: it. {}, obj. {}'.format(it,objective))
+# Initialize conveks
+conveks.initialize()
- # call MMA and update the current point
- conveks.mma.setMoveLimits(lowerBound, upperBound, 0.1)
- conveks.mma.updateCurrentPoint(grad_objective)
+# Initialize the MMA optimizer
+conveks.mma.initialize(initialPoint, lowerBound, upperBound)
+
+# Set some options for MMA
+conveks.mma.option.setNumber('General.Verbosity', 0*4)
+
+# Get iteration count (here it will be 1 - could be different in case of restart)
+it = conveks.mma.getOuterIteration()
+innerit = 0; inneritTot = 0; change = 1.0; kkt_norm = 1.0
+
+# we kkep track of the objective values in all_objectives
+all_objectives = []
+
+# get (copy of) current point
+xFromMMA = conveks.mma.getCurrentPoint()
+
+# compute objective function and constraints
+# as well as their sensitivity with respect to design variables at `x'
+objective,constraints = get_objective(xFromMMA)
+grad_objective,grad_constraints = get_grad()
+
+# send data to onelab server
+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)))
+
+# show current iteration
+print 'iter. inner iter. obj. max(constr.) L2-norm(kkt) change point'
+print '%3i %3i %.4e %.4e %.4e %.4e'%(it-1, innerit, objective, np.max(np.abs(constraints)), kkt_norm, change),
+print xFromMMA[0:5]
+
+while it <= maxIter and change > maxChange \
+ and c.getString('shape/Action') != 'stop':
+
+ # Solve MMA
+ if Optimizer==0:conveks.mma.setMoveLimits(lowerBound, upperBound, 0.1)
+ conveks.mma.updateCurrentPoint(objective,grad_objective,constraints,grad_constraints)
+
+ # 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:
+ # Evaluate the objective at the new point and check if
+ # the approximation is conservative at the current point
+ obj_new,c_new = get_objective(conveks.mma.getCurrentPoint())
+ 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))
+
+ while innerit <= 15 and conserv==0:
+ 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)
+
+ # 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))
+
+ inneritTot += innerit
+
+ # get (copy of) current point
+ xFromMMA = conveks.mma.getCurrentPoint()
+
+ # compute objective function and constraints
+ # as well as their sensitivity with respect to design variables at `x'
+ objective,constraints = get_objective(xFromMMA)
+ grad_objective,grad_constraints = get_grad()
+ all_objectives.append(objective)
+
+ # evaluate the L2-norm of KKT
+ kkt_norm = conveks.mma.getKKTNorm(constraints,grad_objective,grad_constraints)
+
+ # get the change between two successive points
change = conveks.mma.getDesignChange()
+
+ # show the current state
+ print '%3i %3i %.4e %.4e %.4e %.4e'%(it, innerit, objective, np.max(np.abs(constraints)), kkt_norm, change),
+ print xFromMMA[0:5]
+
+ # get the outer iteration count here
it = conveks.mma.getOuterIteration()
+ # send data to the onelab server
+ c.setNumber('Optimization/Results/01Current iteration', value=it, readOnly=1)
+ c.setNumber('Optimization/Results/02Current change', value=change, readOnly=1)
+ c.setNumber('Optimization/Results/03Current KKT norm', value=kkt_norm, readOnly=1)
+ 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)))
+
+# show the problem at the converged state
+print '*'*80
+print 'Optimization has converged after %i iterations.'%(it-1)
+print ' > point :',xFromMMA[0:5]
+print ' > objective : %.4e'%(objective)
+print ' > max(constraints): %.4e'%(np.max(np.abs(constraints)))
+print ' > L2-norm of kkt : %.4e'%(kkt_norm)
+print ' > inner iterations: %i'%(inneritTot)
+print ' > change : %.4e'%(change)
+print '*'*80
+
# This should be called at the end of MMA
conveks.mma.finalize()
--
GitLab