Initialization

This chapter describes the SACOBRA initialization and the options (or hyperparameters) of SACOBRA.

CobraInitializer

The initialization in SACOBRA_Py is done by creating an object of class CobraInitializer and it consists of the following steps:

  • pass in the specification of the COP

  • pass in all the options in SACoptions s_opt

  • (optional, if s_opts.ID.rescale==True) rescale the problem in input space

  • create the initial design, see InitDesigner

  • adjust several elements according to constraint range, see adCon()

  • calculate for each initial design point numViol, the number of violated constraints, and maxViol, the maximum constraint violation. If equality constraints are involved, calculate \(\mu_{init}\), the radius for an artificial feasibility tube around each equality constraint (see EQUoptions) and base the calculation of numViol and maxViol on this artificial feasibility.

  • calculate the so-far best (artificial) feasible point. If no point fulfills (artificial) feasibility, form the set of points with minimum numViol and take from this set the one point with the best objective.

  • set up result dictionary sac_res

  • adjust DRC according to objective range, see adDRC()

class cobraInit.CobraInitializer(x0, fn: object, f_name: str, lower: ~numpy.ndarray, upper: ~numpy.ndarray, is_equ: ~numpy.ndarray, solu=None, s_opts: ~opt.sacOptions.SACoptions = <opt.sacOptions.SACoptions object>)

Initialize SACOBRA optimization:

  • problem formulation: x0, fn, lower, upper, is_equ, solu

  • parameter settings: s_opts via SACoptions

  • create initial design: A, Fres, Gres via InitDesigner

Parameters:

  • x0 (np.ndarray or None) – start point, if given, then its dim has to be the same as lower. If it is/has NaN or None on input, it is replaced by a random point from [lower, upper].

  • fn – function returning (1+nConstraints)-dim vector: [objective to minimize, constraints]

  • f_name – function name

  • lower – lower bound, its dimension defines input space dimension

  • upper – upper bound (same dim as lower)

  • is_equ – boolean vector with dim nConstraints: Which constraints are equality constraints?

  • solu (np.ndarray or None) – (optional, for diagnostics) true solution vector or solution matrix (one solution per row): one or several feasible x that deliver the minimal objective value

  • s_opts (SACoptions) – the options. If not specified, take the default SACoptions object

Objects of class CobraInitializer (usually named cobra below) have the following useful attributes:

  • phase name of the optimization phase [‘init’ | ‘phase1’ | ‘phase2’]

  • rng random number generator

  • rw RescaleWrapper

  • sac_opts object of class SACoptions containing all the options

  • sac_res dictionary with the SACOBRA results from initialization and optimization, see cobra.sac_res in the appendix for details

  • df data frame with diagnostic information from optimization, see cobra.df in the appendix for details

  • df2 data frame with further diagnostic information from optimization, see cobra.df2 in the appendix for details

See Appendix for more details on cobra.sacres, cobra.df, cobra.df2.

adCon()

Adjust several elements according to constraint range.

The following elements of self.sac_res may be changed: ‘fn’, ‘Gres’, ‘Grange’, ‘GrangeEqu’.

adDRC()

Adjust DRC (distance requirement cycle), based on range of Fres

get_fbest()

Return the original objective function value at the best feasible solution.

If no feasible solution was found, return the objective value of the solution with the least maximum violation. (Note: In this case, the objective value may be seemingly ‘better’ than the true optimal objective.)

get_feasible_best()

Return the original objective function value at the best feasible solution.

If no feasible solution was found, return NaN.

get_xbest()
Returns:

best feasible solution in original space

Return type:

np.ndarray

If no feasible solution was found, return the solution with the least maximum violation.

get_xbest_cobra()
Returns:

best feasible solution in COBRA space (maybe rescaled)

Return type:

np.ndarray

is_feasible() → bool
Returns:

True, if a feasible solution has been found

class initDesigner.InitDesigner(x0: ndarray, fn, rng, lower: ndarray, upper: ndarray, is_equ: ndarray, s_opts: SACoptions)

The initial design is a set of P points from input space with dimension d. The problem functions \(f,g,h\) are evaluated at these P points and these evaluated sets form the basis of the later optimization: From the characteristics of the evaluated sets, the classes CobraInitializer and CobraPhaseII deduce decisions about certain adjustments:

CobraPhaseII uses the evaluated sets to train the initial fitness and constraint surrogate models.

In detail, the constructor of InitDesigner does the following: Create the initial design in matrix self.A with shape (P, d) of sample points in (potentially rescaled) input space [lower, upper] \(\subset \mathbb{R}^d\), where P = s_opts.ID.initDesPoints and d = input space dimension. The recipe how to select the sample points is prescribed by the type of initial design s_opts.ID.initDesign.

Apply fn to these points and split the result in objective function (\(f\)) values self.Fres with shape (P,) and constraint function (\(g,h\)) values self.Gres with shape (P,nC), where nC = number of constraints.

Parameters:

  • x0 – For the cases ID.initDesign = "OPTCOBYLA" | "OPTBIASED, take x0 as start point for a short optimization run. In the other ID.initDesign choices, set the last point self.A[-1,:] to x0. x0 is rescaled to [lower, upper] if s_opts.ID.rescale is set.

  • fn – see parameter fn in cobraInit.CobraInitializer

  • rng – RNG (random number generator) from cobraInit.CobraInitializer

  • lower – vector of shape (d,)

  • upper – vector of shape (d,)

  • is_equ – boolean vector with dim nConstraints: Which constraints are equality constraints?

  • s_opts (SACoptions) – the options. Here we use s_opts.cobraSeed and from element IDoptions s_opts.ID the elements initDesign, initBias and initDesPoints.

__call__() → tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]

Return the three results A, Fres and Gres of the initial design

Returns:

tuple (self.A, self.Fres, self.Gres)

Return type:

(np.ndarray, np.ndarray, np.ndarray)

The return tuple contains:

  • self.A: (P,d)-matrix with initial design

  • self.Fres: (P,)-vector with objective evaluated at each initial design point

  • self.Gres: (P,nC)-matrix with constraints evaluated at each initial design point

Types of Initial Design

The initial design creates a matrix self.A with shape (P, d) of P sample points in (potentially rescaled) input space [lower, upper] \(\subset \mathbb{R}^d\), where P = s_opts.ID.initDesPoints and d = input space dimension.

The recipe how to select the sample points is prescribed by s_opts.ID.initDesign:

  • “RANDOM”: uniform random

  • “RAND_REP”: uniform random with reproducible random numbers (both in R and in Python)

  • “LHS”: Latin Hypercube Sampling, see SciPy’s LatinHyperCube

  • “BIASED”: random sample from a normal distribution with mean x0 and standard deviation s_opts.ID.initBias. This is useful if x0 is already in the interesting (e.g. near-feasible or objective-minimizing) region.

  • “OPTCOBYLA”: optimized design: perform a short COBYLA optimization run starting from x0 with the real objective and constraint functions. Extract with the help of FnArchiveFactory from the function evaluations the first P independent points. The idea is that we sample not randomly but in the relevant part of the input space.

  • “OPTBIASED”: establish from a short COBYLA run a new suitable point new_x0 (best feasible or near-feasible). Generate a random sample from a normal distribution with mean new_x0 and standard deviation s_opts.ID.initBias.

class fnArchiveFact.FnArchiveFactory(fn, x0)

Build a wrapper for function fn with input and function value archive.

Usage:

fnArch = FnArchiveFactory(fn, x0)
val = fnArch(x)
#  ...
f_arch = fnArch.getFuncArchive()
Parameters:

  • fn – the function to be wrapped and archived

  • x0 – a typical input vector, needed only to infer the archive sizes

__call__(x)

Call v = fn(x) and store x and function value v in archives.

getFuncArchive()
Returns:

an array where each row is an v from a call to v = fnArch(x)

getSoluArchive()
Returns:

an array where each row is an x from a call to fnArch(x)

Distance Requirement Cycle (DRC)

The Distance Requirement Cycle (DRC) is the vector XI that controls exploration: Each already evaluated infill point is surrounded by a forbidden-sphere of radius XI[c] with c = i mod XI.size (c loops cyclically through XI’s inidices`, that’s where the name cycle comes from). A new infill point is searched under the additional constraint that it has to be a distance XI[c] away from all other already evaluted infill points. The larger XI[c], the more exploration.

XI is set by class SACoptions.

If XI==None, then CobraInitializer will set it, depending on objective range, to short DRC [0.001, 0.0] or long DRC [0.3, 0.05, 0.001, 0.0005, 0.0]. Both vectors contain XI[c] = 0 which enforces exploitation. (If all entries were XI[c] > 0 then a good region would never be exploited further, since the search could never continue in the close vicinity of an already good infill point.)

Options

All paramters (options) for SACOBRA_Py have sensible defaults defined. The user has only to specify those parameters where a setting different from the defaults is desired.

class opt.sacOptions.SACoptions(feval=50, XI=None, skipPhaseI=True, saveIntermediate=False, verbose=1, verboseIter=10, important=True, cobraSeed=42, ID=<opt.idOptions.IDoptions object>, RBF=<opt.rbfOptions.RBFoptions object>, SEQ=<opt.seqOptions.SEQoptions object>, EQU=<opt.equOptions.EQUoptions object>, ISA=<opt.isaOptions.ISAoptions object>, MS=<opt.msOptions.MSoptions object>, RI=<opt.riOptions.RIoptions object>, TR=<opt.trOptions.TRoptions object>)

The collection of all parameters (options) for SACOBRA_Py. Except for some general parameters defined in this class, they are hierarchically organized in nested option classes.

Parameters:

  • feval – number of function evaluations

  • XI – Distance-Requirement-Cycle (DRC) that controls exploration: Each infill point has a forbidden-sphere of radius XI[c] around it. c loops cyclically through XI’s inidices. If XI==None, then CobraInitializer will set it, depending on objective range, to short DRC [0.001, 0.0] or long DRC [0.3, 0.05, 0.001, 0.0005, 0.0].

  • skipPhaseI – whether to skip SACOBRA_Py phase I or not

  • saveIntermediate – whether to save intermediate results (including surrogates) or not. If True, save to f'results/cobra-{f_name}-{SEQ.optimizer}-{cobraSeed}.pkl'

  • verbose – verbosity level: 0: print nothing. 1: print only important messages. 2: print every message

  • verboseIter – an integer value, after how many iterations to print summarized results.

  • important – controls the importance level for some verboseprint’s in updateInfoAndCounters

  • cobraSeed – the seed for RNGs. SACOBRA_Py guarantees the same results for the same seed

  • ID (IDoptions) – nested options for initial design

  • RBF (RBFoptions) – nested options for radial basis functions

  • SEQ (SEQoptions) – nested options for sequential optimizer

  • EQU (EQUoptions) – nested options for equality constraints

  • ISA (ISAoptions) – nested Internal SACOBRA options

  • MS (MSoptions) – nested options for model selection (TODO)

  • RI – nested options for repair infeasible (TODO)

  • TR (TRoptions) – nested options for trust region (TODO)

class opt.idOptions.IDoptions(initDesign='RANDOM', initDesPoints: int | None = None, initDesOptP: int | None = None, initBias=0.005, rescale=True, newLower=-1, newUpper=1)

Options for the initial design (\(d\) = input dimension of problem).

Parameters:

  • initDesign – options: “RANDOM”, “RAND_R”, “RAND_REP”, “LHS”, … (see Types of initial design)

  • initDesPoints – number of initial design points. If None, cobraInit will set it to \(d+1\) if RBF.degree=1 or to \((d+1)(d+2)/2\) if RBF.degree=2

  • initDesOptP – if None, cobraInit will set it to initDesPoints

  • initBias – standard deviation for normal distribution in cases initDesign = "BIASED" | "OPTBIASED"

  • rescale – if True, rescale input space from [lower, upper] to \([\) newLower, newUpper \(]^d\)

  • newLower – common new lower bound for each of the \(d\) input dimensions

  • newUpper – common new upper bound for each of the \(d\) input dimensions

class opt.rbfOptions.RBFoptions(kernel='cubic', degree=2, rho=0.0, rhoDec=2.0, rhoGrow=0, width=None, widthFactor=1.0, widthRule=W_RULE.ONE, interpolator='scipy', test_pmat=False)

Options for the RBF surrogate models

Parameters:

  • kernel – RBF kernel type, see below

  • degree – degree of polynomial tail for RBF kernel. See SciPy’s RBFInterpolator for details.

  • rho – Smoothing parameter \(\rho\). If \(\rho=0\), use interpolating RBFs: The surrogate model surface passes exactly through the points. If \(\rho>0\), use approximating RBFs (spline-like). The larger \(\rho\), the smoother the surrogate model.

  • rhoDec – exponential decay factor for \(\rho\)

  • rhoGrow – every rhoGrow (e.g. 100) iterations, re-enlarge \(\rho\). If 0, then re-enlarge never

  • width –

    only relevant for the scale-variant kernel types. width = \(\sigma\) translates to shape parameter epsilon = \(\epsilon = 1/\sqrt{2\sigma}\) in SciPy’s RBFInterpolator.

  • widthRule – only relevant for the scale-variant kernels: If width=None, calculate the appropriate width from the data by heuristic rule W_RULE.ONE or W_RULE.THREE

  • widthFactor – only for scale-variant kernels. Additional constant factor applied to each width \(\sigma\)

  • interpolator –

    ”scipy” or “sacobra”, which interpolation method to use. In case of “scipy”, use SciPy’s RBFInterpolator (faster and simpler in code). In case of “sacobra”, use an object of class RBFsacob, which is SACOBRA’s own implementation of RBF models (allows with degree=1.5 the option equivalent to squares=T in SACOBRA R, which means only pure squares in the polynomial tail).

  • test_pmat – only for testing the RBFsacob implementation (degree=1, 1.5)

Parameter kernel specifies the kernel type and is one out of "cubic", "quintic", "thin_plate_spline", "gaussian", "multiquadric".

  • scale-invariant kernel types: "cubic", "quintic", "thin_plate_spline"

  • scale-variant kernel types: "gaussian", "multiquadric"

class opt.seqOptions.SEQoptions(optimizer='COBYLA', feMax=1000, tol=1e-06, conTol=0.0, penaF=[3.0, 1.7, 300000.0], sigmaD=[3.0, 2.0, 100], epsilonInit=None, epsilonMax=None, finalEpsXiZero=True, Tfeas=None, Tinfeas=None, trueFuncForSurrogates=False)

Options for the sequential optimization

Parameters:

  • optimizer – string defining the optimization method for SACOBRA phase I and II. One out of [“COBYLA”,”ISRESN”] (see SeqOptimizer)

  • feMax – maximum number of function evaluations on the surrogate model

  • tol – convergence tolerance for sequential optimizer

  • conTol – constraint violation tolerance

  • penaF – (TODO)

  • sigmaD – (TODO)

  • epsilonInit – initial value for EPS. EPS is a constant added to each constraint to maintain a certain margin to the boundary. If None, then epsilonInit is set to \(0.005 \ell\), where \(\ell\) is the length of smallest side of search space

  • epsilonMax – maximum value for EPS. If None, then epsilonMax is set to \(2*0.005 \ell\).

  • finalEpsXiZero – if True, set in final iteration EPS and XI to zero for best exploitation

  • Tfeas – threshold for count of feasible iterations in a row. If None, CobraInitializer will set it to \(floor(2\sqrt{d})\).

  • Tinfeas – threshold for count of infeasible iterations in a row. If None, CobraInitializer will set it to \(floor(2\sqrt{d})\).

  • trueFuncForSurrogates – if True, use the true (constraint & fitness) functions instead of surrogates (only for debug analysis)

class opt.equOptions.EQUoptions(active=True, initType='TAV', muType='expFunc', muDec=1.5, muFinal=1e-07, mu4inequality=False, muGrow=0, refine=True, refineMaxit=1000, refineAlgo='BFGS_0', refinePrint=False)

Equality handling options

Parameters:

  • active – if set to TRUE, the equality-handling (EH) technique is activated. The EH technique transforms each equality constraint \(h(\mathbf{x})=0\) into two inequality constraints \(h(\mathbf{x})-\mu <0\) and \(-h(\mathbf{x})-\mu<0\) with an adaptive (normally decaying) margin \(\mu\).

  • initType – the equality margin \(\mu\) is initialized with one of these choices: ["TAV" | "TMV" | "EMV" | "useGrange"]

  • muType – type of function used to shrink margin \(\mu\) during the optimization process. One out of ["SAexpFunc" | "expFunc" | "funcDim" | "funcSDim" | "Zhang" | "CONS"], see modifyMu in equHandling.py

  • muDec – decay factor for margin \(\mu\), see modifyMu

  • muFinal – lower bound for margin \(\mu\). muFinal should be set to a small but non-zero value (larger than machine accuracy).

  • muGrow – every muGrow (e.g. 100) iterations, re-enlarge the \(\mu\)-band. If 0, then re-enlarge never

  • mu4inequality – use the artificial feasibility band also for inequalities (experimental)

  • refine – enables the refine step for equality handling

  • refineMaxit – maximum number of iterations used in the refine step. Note that the refine step runs on the surrogate models and does not impose any extra real function evaluations

  • refineAlgo – optimizer for refine step [“BFGS_0” | “BFGS_1” | “COBYQA” | “COBYLA”]

  • refinePrint – whether to print “cg-values (before,after,true)” after each refine step

class opt.riOptions.RIoptions(repairInfeas=True, eps1=0.0001, eps2=0.0001, q=3.0, mmax=1000, gradEps=None, repairMargin=0.01, repairOnlyFresBetter=False, marFres=0.0, trueFuncForSurrogates=False, checkIt=False)

Options for RI2 (method to repair infeasible solutions)

The infeasibility of a solution is its maximum constraint violation (0 for a feasible solution).

A solution \(x\) is \(\epsilon\)-feasible for constraint function \(g\) if \(g(x) + \epsilon < 0\).

Parameters:

  • repairInfeas – switch to turn repair infeasible on or off

  • eps1 – include all constraints not eps1-feasible into the repair mechanism

  • eps2 – selects the solution with the shortest shift among all random realizations which are eps2-feasible

  • q – draw coefficients \(\alpha_k\) from uniform distribution \(U[0, q]\)

  • mmax – draw mmax random realizations

  • gradEps – stepsize for numerical gradient calculation. If None, set it to 1/1000 of smallest search space box side

  • repairMargin – repair only solutions whose infeasibility is less than this margin

  • repairOnlyFresBetter – if True, then repair only iterates with fitness < so-far-best-fitness + marFres

  • marFres – only relevant if repairOnlyFresBetter==True

  • trueFuncForSurrogates – use true constraint functions instead of constraint surrogates

  • checkIt – if true, print debug and check info

class opt.isaOptions.ISAoptions(isa_ver=1, RS=True, RStype=RSTYPE.CONSTANT, RSauto=False, RSmax=0.3, RSmin=0.05, RS_Cs=10, RS_rep=False, RS_verb=False, aDRC=True, aFF=True, aCF=True, TFRange=100000.0, TGR=1000.0, adaptivePLOG=False, onlinePLOG=O_LOGIC.NONE, onlineFreqPLOG=10, pEffectInit=0, pEffNpts=3, pEff_DBG=False)

Internal self-adjusting (ISA) options for SACOBRA.

Defaults for isa_ver=1 (SACOBRA settings with fewer online adjustments). Other choices are

  • isa_ver=0 (ISAoptions0, COBRA-R settings = turn off SACOBRA), and

  • isa_ver=2 (ISAoptions2, SACOBRA settings with fewer parameters and more online adjustments).

The choice is made by choosing the appropriate class object ISAoptions, ISAoptions0 or ISAoptions2 for cobra.sac_opts.ISA.

Parameters:

  • isa_ver – ISA version number (was DOSAC in R)

  • RS – flag to enable/disable random start algorithm

  • RStype – type of function to calculate probability to start the internal optimizer with a random start point. One out of [RSTYPE.CONSTANT, RSTYPE.SIGMOID]

  • RSauto –

  • RSmax – maximum probability of a random start

  • RSmin – minimum probability of a random start

  • RS_Cs – if RS_Cs iterations w/o progress, do a random start

  • RS_rep – if True, generate reproducible random numbers with my_rng2 (R and Python)

  • RS_verb – if True, be verbose in RandomStarter.random_start

  • aDRC – flag for automatic DRC adjustment

  • aFF – flag for automatic objective function transformation

  • aCF – flag for automatic constraint function transformation

  • TFRange – threshold, if the min-max-range of Fres[0:idp] is larger than TFRange and if onlinePLOG=NONE, then apply plog to Fres (objective function values)

  • TGR – threshold: If GRatio > TGR, then apply automatic constraint function transformation. GRatio is the ratio “largest GR / smallest GR” where GR is the min-max range of a specific constraint. If TGR < 1, then the transformation is always performed.

  • adaptivePLOG – (experimental) flag for objective function transformation with plog, where the parameter pShift is adapted during iterations

  • onlinePLOG – three-valued logic for online decision-making: O_LOGIC.NONE (no online, only fixed initial decision whether to use plog); O_LOGIC.XNEW (online plog decision according to pEffect, where pEffect-calculation is based on new infill point xNew); O_LOGIC.MIDPTS (online plog decision according to pEffect, where pEffect-calculation is based on midpoints)

  • onlineFreqPLOG – after how many iterations the surrogates for online plog check are calculated again (only relevant for case O_LOGIC.XNEW)

  • pEffectInit – the initial value for pEffect, needed for first pass through cobraPhaseII while loop in case O_LOGIC.XNEW. Not needed in cases O_LOGIC.NONE or O_LOGIC.MIDPTS

  • pEffNpts – only for case O_LOGIC.MIDPTS: the number p of initial design points to take from initial design matrix A: Form p*(p-1)/2 pairs and calculate the midpoints for all those pairs

  • pEff_DBG – temporary debug: replace xNew with a single midpoint (needed only by testCOP.py, which activates a certain branch in cobraPhaseII.py)

class opt.isaOptions.ISAoptions0(isa_ver=0, RS=False, RStype=RSTYPE.CONSTANT, RSauto=False, RSmax=0.3, RSmin=0.05, RS_Cs=10, RS_rep=False, RS_verb=False, aDRC=False, aFF=False, aCF=False, TFRange=inf, TGR=inf, adaptivePLOG=False, onlinePLOG=O_LOGIC.NONE, onlineFreqPLOG=10, pEffectInit=0, pEffNpts=3, pEff_DBG=False)

Internal self-adjusting options for isa_ver=0 (COBRA-R settings = turn off SACOBRA).

Differs only in default settings. The parameters and their meaning are the same as in ISAoptions.

class opt.isaOptions.ISAoptions2(isa_ver=2, RS=True, RStype=RSTYPE.CONSTANT, RSauto=True, RSmax=0.3, RSmin=0.05, RS_Cs=10, RS_rep=False, RS_verb=False, aDRC=True, aFF=True, aCF=True, TFRange=-1, TGR=-1, adaptivePLOG=False, onlinePLOG=O_LOGIC.NONE, onlineFreqPLOG=10, pEffectInit=3, pEffNpts=3, pEff_DBG=False)

Internal self-adjusting options for isa_ver=2 (SACOBRA with fewer parameters and more online adjustments).

Differs only in default settings. The parameters and their meaning are the same as in ISAoptions.