Simulation Optimization Methods - Part 2

Learning Objectives

• To be able discuss common behaviors and properties of simulation optimization solvers.
• To be able to describe the iterative nature of solvers.
• To be able to represent adaptable solver behavior
• To be able to understand how solutions are managed in solvers.
• To be able to discuss and use specific solvers

Characterizing Simulation Optimization Algorithms

• The KSL ksl.simopt.solvers package provides a general framework for designing and implementing simulation optimization algorithms.
  
• We call these algorithms, solvers.
  
• The solver package provides the Solver abstract base class that can be specialized to implement different algorithms.

The Base Solver Class

abstract class Solver(
    val problemDefinition: ProblemDefinition,
    evaluator: EvaluatorIfc,
    maximumIterations: Int,
    var replicationsPerEvaluation: ReplicationPerEvaluationIfc,
    name: String? = null
) : IdentityIfc by Identity(name), Comparator<Solution>, SolverEmitterIfc by SolverEmitter() {

• A solver must have a reference to the problem and an evaluator.
• Solvers are conceptualized as iterative processes. Each iteration of a solver may evaluate solutions and move towards a “best” solution.
  • Thus, solvers require a limit on the maximum number of iterations.
• Because of the stochastic nature of the evaluation, the solver may request one or more replications for its evaluation requests.
  • The number of replications may change dynamically, and thus the user needs to supply a function to determine the number of replications per evaluation.

Stochastic Solvers

The StochasticSolver inherits from the Solver class and requires that a source of random streams be provided. It also allows stream control.

abstract class StochasticSolver(
    problemDefinition: ProblemDefinition,
    evaluator: EvaluatorIfc,
    maxIterations: Int,
    replicationsPerEvaluation: ReplicationPerEvaluationIfc,
    streamNum: Int = 0,
    val streamProvider: RNStreamProviderIfc = KSLRandom.DefaultRNStreamProvider,
    name: String? = null
) : Solver(problemDefinition, evaluator, maxIterations, replicationsPerEvaluation, name), RNStreamControlIfc {

Solvers as Iterative Processes

A solver is an iterative algorithm that searches for the optimal solution for a defined problem. In the abstract Solver class, the algorithm is conceptualized as having a main iterative loop. The basic algorithm is structured conceptually as an iterative loop:

initializeIterations()
while(hasNextIteration()){
  mainIteration()
  iterationCounter++
  checkStoppingCondition()
}
endIterations()

Required Abstract Methods

The Solver class has three abstract protected functions that must be implemented.

• protected abstract fun startingPoint(): InputMap - The purpose of this function is to determine the starting point for the solver’s process. In the default implementation, this function is called from the initializeIterations() function.
• protected abstract fun nextPoint(): InputMap - The purpose of this function is to determine the next point to be evaluated for the solver’s process. Specific implementations will call this function within their implementations of the mainIteration() function.
• protected abstract fun mainIteration() - The purpose of this function is to contain the logic that iteratively executes until the maximum number of iterations is reached or until the stopping criteria is met.

Default Iteration Functionality

The Solver class has protected open functions that can be overridden with more specific behavior:

• protected open fun initializeIterations() - This function is called prior to the running of the main loop of iterations.
  • The base implementation determines a starting point, requests an evaluation, and sets the initial current solution.
• protected open fun beforeMainIteration() - An empty hook function to permit logic to be inserted prior to the call to mainIteration().
• protected open fun afterMainIteration() - An empty hook function to permit logic to be inserted after the call to mainIteration().
• protected open fun mainIterationsEnded() - An empty hook function to permit logic to be inserted after the main iterative loop completes.

Default Stopping Functionality

• protected open fun isStoppingCriteriaSatisfied(): Boolean - This function is called from checkStoppingCondition() within the iterative loop.
• Its purpose is to implement a stopping criteria that is based on the quality of the solution.
• The user has the option to specify a function via the solutionQualityEvaluator property that will be used from within this function.
• If no solution quality checker is supplied or this function is not overridden, then this function should report false.
• In that case, the main loop will only stop when the maximum number of iterations has been reached.

Generating Neighbors

• protected open fun generateNeighbor(currentPoint: InputMap,rnStream: RNStreamIfc): InputMap
• A basic function that can be used by sub-classes to generate a point that is a neighbor to the current point. The process can be random.
• The user has the option of either overriding this function or supplying an instance of the GenerateNeighborIfc interface via the neighborGenerator property.
• If a neighborGenerator is not supplied, the base algorithm, randomly selects one of the coordinates from the current point and then randomly generates an input range feasible value for the coordinate.
  • An input range feasible value only checks the bounds on the decision variable and does not utilize linear or functional constraints.

Stochastic Hill Climbing Algorithm

• This algorithm searches for an optimal solution by iteratively evaluating and moving to neighborhood solutions, as determined by the evaluator and the problem definition.
• It uses a stochastic approach to explore the solution space by leveraging a random number stream.
• Randomly generated points are evaluated and if the new point has a better solution than the current solution, the proposed point becomes the current solution.
• This process continues until the maximum number of iterations is met.

Stochastic Hill Climbing Main Iteration Function

    override fun mainIteration() {
        // generate a random neighbor of the current solution
        val nextPoint = nextPoint()
        // evaluate the solution
        val nextSolution = requestEvaluation(nextPoint)
        if (compare(nextSolution, currentSolution) < 0) {
            currentSolution = nextSolution
        }
        // capture the last solution
        solutionChecker.captureSolution(currentSolution)
    }

Stochastic Hill Climbing Starting and Next Points

• The implementations of startingPoint() and nextPoint() both inherited from StochasticSolver are as follows:

    override fun startingPoint(): InputMap {
        return startingPointGenerator?.startingPoint(problemDefinition) ?: problemDefinition.startingPoint(rnStream)
    }
    
    override fun nextPoint(): InputMap {
        return generateNeighbor(currentPoint, rnStream)
    }

• nextPoint() leverages the availability of the supplied generateNeighbor() function.
• The implementation of the startingPoint() function will use an instance of the StartingPointIfc interface via the startingPointGenerator property or rely on the problem definition to randomly generate an input feasible point.
  • An input feasible point is a point in the search space that is feasible with respect to the input ranges, linear constraints, and the functional constraints.

Basic Solver Requirements

An implementation of a solver must at the very least specify:

• how to determine its starting point (startingPoint()),
• how to determine its next point for evaluation (nextPoint()), and
• how to search its the solution space (mainIteration()).

New solvers can be constructed by sub-classing from Solver or StochasticSolver or by modifying the basic behavior by plugging in different functions.

Solver Plug-in Behavior

Functional interfaces provide the ability to specify plug-in behavior by supplying a different implementation for one of these interfaces.

• ReplicationPerEvaluationIfc interface - This functional interface can be used by solvers to determined the number of replications to request when asking the simulation oracle for evaluations.
• GenerateNeighborIfc interface - This functional interface can be used by solvers to specify how to generate a neighbor relative to the supplied input.
• NeighborSelectorIfc interface - This functional interface defines functions for selecting a neighbor (point) from a defined neighborhood represented as as set containing inputs.
• SolutionQualityStoppingCriteriaIfc interface - This functional interface represents functions that define when to stop the iterative process based on the quality of the solution.

Replications Per Evaluation

fun interface ReplicationPerEvaluationIfc {

    fun numReplicationsPerEvaluation(solver: Solver) : Int
}

• This functional interface can be used by solvers to determined the number of replications to request when asking the simulation oracle for evaluations.
• This permits strategies that start with less replications during the initial search and then increase replications when exploring higher quality solutions.
  
• There are two concrete implementations: the FixedReplicationsPerEvaluation and FixedGrowthRateReplicationSchedule classes.
  • The fixed schedule causes a constant number of replications to be used for every evaluation request.
  • The fixed growth rate replication schedule models an increasing number of replications based on a growth rate factor that is a function of the number of iterations.

Replications Per Evaluation

Search Neighborhoods

• GenerateNeighborIfc interface - This functional interface can be used by solvers to specify how to generate a neighbor relative to the supplied input.
  
• The ensureFeasibility parameter should indicate if the generation method should ensure that the returned input map is feasible with respect the constraints of the problem.

fun interface GenerateNeighborIfc {

    fun generateNeighbor(inputMap: InputMap, solver: Solver, ensureFeasibility: Boolean) : InputMap

}

Specifying Your Own Neighborhoods

• The current generateNeighbor() function randomly select one of the coordinates (inputs) and then randomly generates an input-range feasible value for the selected input.
  • The randomly generated point may be infeasible with respect to deterministic constraints.
• To specify your own neighborhood generation procedure, either
  • override the generateNeighbor(currentPoint: InputMap,rnStream: RNStreamIfc) function, or
  • implement the GenerateNeighborIfc interface and set the neighborGenerator property with the instance

Selecting a Point from a Neighborhood

• NeighborSelectorIfc interface - This functional interface defines functions for selecting a neighbor (point) from a defined neighborhood represented as as set containing inputs. The supplied solver could be used to allow memory during the selection process.

fun interface NeighborSelectorIfc {

    fun select(neighborhood: Set<InputMap>, solver: Solver) : InputMap

}

Specifying Stopping Criteria

• SolutionQualityStoppingCriteriaIfc interface - This functional interface represents functions that define when to stop the iterative process based on the quality of the solution.

• The solutionQualityEvaluator property will be used when the stopping criteria is evaluated within the main loop.

• The solutionQualityEvaluator property provides the user with an approach to overriding the base behavior with different behavior rather than through sub-classing.

fun interface SolutionQualityEvaluatorIfc {

    fun isStoppingCriteriaReached(solver: Solver): Boolean

}

Other Solver Functionality

Key Properties

• currentPoint: Represents the current point (input settings) of the solver during its iterative process.
  • This value is derived from the inputMap of the currently active solution, which reflects the solver’s progress in finding optimal or improved solutions.
• currentSolution: The current (or last) solution that was accepted as a possible solution to recommend for the solver. It is the responsibility of the subclass to determine the current solution.
• initialSolution: The initial starting solution for the algorithm. It is the responsibility of the subclass of Solver to initialize the initial solution.
• startingPoint: An initial starting point for the solver. If supplied, this point will be used instead of the returned value of the startingPoint() function.
  • The default is null, which indicates that the function should be called to obtain the initial starting point.
  • The starting point must be a valid point in the input space.
  • It must also be input-feasible.
• previousPoint and previousSolution: The previously evaluated point and the previously resulting solution.

Other Properties

• iterationCounter - The iterationCounter property counts the number of iterations performed by a solver. Since some solvers may used additional stopping criteria, this property can be used to assess the total number of iterations performed.
• numOracleCalls - Each time a solver requests evaluations from the oracle, this property is incremented. It can be used to indicate how may simulation evaluations were performed. Since simulation evaluations can be computationally expensive, this property can help understand how efficient a solver is with respect to simulation evaluations.
• numReplicationsRequested - Similar to the numOracleCalls property, this property totals the number of replications requested during the solver’s execution. This property can help understand how efficient a solver is in terms of its sampling.

Underlying Iterative Process

Understanding Solutions

• currentSolution: The current (or last) solution that was accepted as a possible solution to recommend for the solver.
  • The current solution is initialized to a bad solution as defined by the problem.

    var currentSolution: Solution = problemDefinition.badSolution()
        protected set(value) {
            // save the previous solution
            previousSolution = field
            // update the current solution
            field = value
            myBestSolutions.add(value)
        }

• The myBestSolutions protected property is used to capture the newly assigned solution. The protected myBestSolutions property is an instance of the Solutions class.
• The bestSolutions property provides access to the solutions via the SolutionsIfc interface.

Capturing Solutions

Solution Checkers

• The SolutionChecker class is designed to hold and test if its contained solutions are all the same.
• The purpose of the SolutionEqualityIfc interface is to provide a mechanism to check for the equality between two solutions

Checking for Equality

• ConfidenceIntervalEquality - This class checks for equality between solutions based whether the confidence interval on the difference contains the indifference zone parameter. This is based on an approximate confidence interval on the difference between the unpenalized objective function values.
• PenalizedObjectiveFunctionEquality - This class checks if the penalized objective function values are within some specified precision. The uncertainty of the values is not taken into consideration.
• InputEquality This class indicates that two solutions are equal based only on the values of their inputs. If the solutions do not have the same inputs, then the penalized objective function is used to determine the ordering. Thus, two solutions are considered the same if they have the same input values, regardless of the value of the objective functions.
• InputsAndConfidenceIntervalEquality - This class checks for equality between solutions based whether the confidence interval on the difference contains the indifference zone parameter and whether the input variable values are the same.

The SolutionChecker Class

• The SolutionChecker class uses an instance of the SolutionEqualityIfc interface to check if the elements that it holds all test as equal.
• The constructor requires a threshold value that represents the number of solutions to check.
• When the solution checker contains the specified number of solutions, the user can call the checkSolutions() function, which will return true if the all the contained solutions test as equal. The specified threshold limit also specifies the capacity.
• The checker will only hold the specified number of solutions.

class SolutionChecker(
    var equalityChecker: SolutionEqualityIfc = InputEquality,
    noImproveThreshold: Int = defaultNoImproveThreshold
)

• An instance of the SolutionChecker class can be used as part of the logic for checking if the stopping criteria for the algorithm indicates if the algorithm should stop.

Simulated Annealing

• Simulated annealing exploits the MCMC sampling process to find a mode of density function, \(f(x)\). This mode represents a point where \(f(x)\) is maximized.
• The process produces a sequence of density functions \(f_T(x) \propto [f(x)]^{1/T}\), where \(T\) denotes the temperature of the distribution.
  
• The sampling process produces a Markov chain such that the generated values \(x_k\) are from the distributions \(f_{T_k}\) for decreasing temperatures \(T_1, T_2, \dots\).
  • As the temperatures decrease, the distributions converge to a highly peaked distribution, where the peak is at the maximum of the distribution, \(f\).
  • The point causing the peak will also be the solution to the optimization problem.
  • A key factor in this process is the cooling schedule that determines the temperature sequence.

Overview of Simulated Annealing

Pseudo-code for the overall algorithm is as follows:

1. Initialize the solution and the temperature
2. While not converged
  2.1. Generate a possible solution based on some probabilistic process
  2.2. Calculate the change in the objective function
  2.3. If an improvement or acceptance probability is met
    2.3.1 accept the new solution
  2.4 Reduce the temparture according the cooling schedule
3. Return the best solution found

Implementing Simulated Annealing

    override fun mainIteration() {
        // generate a random neighbor of the current solution
        val nextPoint = nextPoint()
        // evaluate the point to get the next solution
        val nextSolution = requestEvaluation(nextPoint)
        // calculate the cost difference
        costDifference = nextSolution.penalizedObjFncValue - currentSolution.penalizedObjFncValue
        // if the cost difference is negative, it is automatically an acceptance
        if (costDifference < 0.0) {
            // no need to generate a random variate
            currentSolution = nextSolution
            lastAcceptanceProbability = 1.0
        } else {
            // decide whether an increased cost should be accepted
            val u = rnStream.randU01()
            lastAcceptanceProbability = acceptanceProbability(costDifference, currentTemperature)
            if (u < lastAcceptanceProbability) {
                currentSolution = nextSolution
            }
        }
        // update the current temperature based on the cooling schedule
        currentTemperature = coolingSchedule.nextTemperature(iterationCounter)
        // capture the last solution
        solutionChecker.captureSolution(currentSolution)
    }

Simulated Annealing Acceptance Probability

The acceptance probability criteria is based on the Boltzman distribution and is a function of the objective function difference and the current value of the temperature.

    fun acceptanceProbability(costDifference: Double, temperature: Double) : Double {
        require(temperature > 0.0) { "The temperature must be positive" }
        return if (costDifference <= 0.0) {
            1.0
        } else {
            exp(-costDifference / temperature)
        }
    }

Cooling Schedule

• The current temperature is computed based on a user definable function to implement the cooling schedule.
  
• The default implementation uses an exponential cooling schedule.
  
• The ksl.solvers.algorithms package also provides implementations for linear and logarithmic cooling schedules.

class ExponentialCoolingSchedule(
    initialTemperature: Double,
    val coolingRate: Double = defaultCoolingRate
) : CoolingSchedule (initialTemperature)  {
    init {
        require(coolingRate > 0.0) { "Cooling rate must be positive" }
        require(coolingRate < 1.0) { "Cooling rate must be less than 1.0" }
    }

    override fun nextTemperature(iteration: Int): Double {
        require(iteration > 0) { "The iteration must be positive. It was $iteration." }
        return initialTemperature * coolingRate.pow(iteration)
    }
}

Stopping Criteria

• The implementation of the SimulatedAnnealing class uses an instance of the previously discussed SolutionChecker class as part of the stopping criteria to check if the last 5 (default) solutions are the same or if the temperature has become too small.
• The user can specify number of solutions to check.

    override fun isStoppingCriteriaSatisfied(): Boolean {
        return solutionQualityEvaluator?.isStoppingCriteriaReached(this) ?:
        checkTemperature() || solutionChecker.checkSolutions()
    }
    
    private fun checkTemperature() : Boolean {
        return currentTemperature < stoppingTemperature
    }

Cross Entropy Solver

• The cross-entropy method is a stochastic optimization procedure that can be applied to both deterministic and stochastic optimization problems.
• The basic procedure can be characterized as a population based approach because it uses a sample of the search space to define an elite population.
• Cross-entropy uses rare-event simulation techniques to cause the parameters of a sampling distribution to converge to the desired solution of the related optimization problem.

Embedded Stochastic Program

• The stochastic program mentioned in step 2.5 of the algorithm is essentially the estimation of the parameters of \(f(\cdot;v_{k})\) via the maximum likelihood method based on the elite sample at each iteration.
• The elite sample consists of those \(\vec{X}_i\) such that \(G(\vec{X}_i) \ge \hat{\gamma}_k\).
• Implementations recommend that the procedure be stopped if for some \(k \ge d\), say \(d=5\), when \(\hat{v}_k = \hat{v}_{k-1} = \cdots = \hat{v}_{k-d}\).
• The choice of the density function \(f(\cdot;v_{k})\) is an important aspect of the algorithm because the parameters converge to the optimal solution of the associated optimization problem.
  • Actually, similar to simulated annealing, the sequence of density functions \(f(\cdot;\hat{v}_{1}),f(\cdot;\hat{v}_{2}), \dots\) converges to a degenerate distribution assigning all probability to a single point \(x_k\), for which the objective function value is greater than or equal to \(\hat{\gamma}_k\). This point represents the solution to optimization problem.
• Because of the ease of solving for the MLE for the normal distribution it is the recommended default distribution.

KSL Cross Entropy Solver

• The KSL implementation of the cross-entropy solver follows this general algorithm and builds from the StochasticSolver class.
• You must specify an instance of the CESampler class or accept the default of a sampler based on the multivariate normal distribution.
• The stopping criteria is based on providing an instance of the InputsAndConfidenceIntervalEquality class, which implements the SolutionEqualityIfc interface.
  • This will test if \(\hat{v}_k = \hat{v}_{k-1} = \cdots = \hat{v}_{k-d}\) is true based on equal values of inputs and a comparison of the solutions based on confidence interval comparisons with the best solution.
• The implementation uses the instance of the CESampler class to generate the points \(\vec{X}_1, \vec{X}_2, \dots, \vec{X}_N\) via the function call ceSampler.sample(sampleSize()) within the mainIteration() function.
• The function sampleSize() computes the size, \(N\), of the cross-entropy population at each iteration.
• Then, the solver requests the evaluation of those points to produce the estimated performance of the objective function, \(G(\vec{X}_1), G(\vec{X}_2),\dots, G(\vec{X}_N)\).

KSL Cross Entropy Solver Constructor

class CrossEntropySolver @JvmOverloads constructor(
    problemDefinition: ProblemDefinition,
    evaluator: EvaluatorIfc,
    streamNum: Int = 0,
    streamProvider: RNStreamProviderIfc = RNStreamProvider(),
    ceSampler: CESampler = CENormalSampler(problemDefinition),
    maxIterations: Int = ceDefaultMaxIterations,
    replicationsPerEvaluation: ReplicationPerEvaluationIfc,
    solutionEqualityChecker: SolutionEqualityIfc = InputsAndConfidenceIntervalEquality(),
    name: String? = null
) : StochasticSolver(problemDefinition,
    evaluator, maxIterations,
    replicationsPerEvaluation, streamNum,
    streamProvider, name
)

Cross Entropy Main Iteration Function

    override fun mainIteration() {
        // Generate the cross-entropy sample population
        val points = ceSampler.sample(sampleSize())
        // Convert the points to be able to evaluate.
        val inputs = problemDefinition.convertPointsToInputs(points)
        // request evaluations for solutions
        val evaluations = requestEvaluations(inputs)
        val results = evaluations.values.toList()
        if (results.isEmpty()) {
            // Returning will cause no updating on this iteration.
            // New points will be generated for another try on the next iteration.
            return
        }
        // At least one result, so proceed with processing.
        // Process the results to find the elites. This should fill myEliteSolutions.
        myEliteSolutions = findEliteSolutions(results)
        // convert elite solutions to points
        val elitePoints = extractSolutionInputPoints(myEliteSolutions)
        // update the sampler's parameters
        ceSampler.updateParameters(elitePoints)
        // specify the current solution
        currentSolution = myEliteSolutions.first()
        // capture the last solution
        solutionChecker.captureSolution(currentSolution)
    }

Meta Parameters for Cross Entropy

There are a number of key issues that must be implemented for the cross-entropy method.

• How large should the population (\(N\)) be?
• What should the elite percentage \(\rho\) be? Alternatively, what should be the size of the elite sample?
• What parametric family of distributions should be used for the sampling distribution?
  • How to initialize the parameters of the sampling distribution?
  • How to update the parameters of the sampling distribution?
  • What should the value of the smoothing parameter \(\alpha\) be?

Sample Size Determination

• The KSL provides more flexibility with respect to determining the population size. The base implementation will use a supplied sample size function or use the value of the ceSampleSize property.

    fun sampleSize(): Int {
        return sampleSizeFn?.sampleSize(this) ?: ceSampleSize
    }

• The value of the ceSampleSize property is based on a user changeable default value or a recommended sample size based on the size of a sample needed to adequately estimate the quantile of a distribution.

    var ceSampleSize: Int = minOf(recommendCESampleSize(), defaultMaxCESampleSize)
        set(value) {
            require(value >= defaultMinCESampleSize) { "The CE sample size must be >= $defaultMinCESampleSize" }
            field = minOf(value, defaultMaxCESampleSize)
        }

Population Size

• The size of the population, \(N\), should be related to the dimension of the problem.
• If we let \(K\) be the number of parameters that need to be estimated for the sampling distribution, then standard implementations recommend \(N=cK\), where \(c\) is a number between 1 and 10.
  • For example, if the sampling distribution family is the normal distribution, we must estimate two parameters for each input variable, that is \((\mu_i, \sigma_i)\) for \(i=1,2,\dots, |I|\).
  • Thus, \(K=2*|I|\).

Size of Elite Sample

• The KSL cross-entropy implementation permits customization of the size of the elite sample.
• In some situations, you may want to vary the size of the elite sample as the algorithm progresses to improve estimation of the parameters of the sampling distribution.
• The KSL implementation provides the eliteSize() function for this purpose.

    fun eliteSize(): Int {
        return eliteSizeFn?.eliteSize(this) ?: maxOf(ceil(elitePct * ceSampleSize).toInt(),
            defaultMinEliteSize)
    }

• Note that a user can supply an instance of the EliteSizeIfc functional interface, which simply promises to compute the size of the elite sample.

fun interface EliteSizeIfc {

    fun eliteSize(ceSolver: CrossEntropySolver): Int

}

• The implementation of the eliteSize() function will use a supplied function specified by the eliteSizeFn property or determine the size based on the specified elite percentage.

The Cross Entropy Sampler

The questions associated with the sampling distribution and the updating of its parameters should be addressed by the choice of the implementation of the CESampler class chosen for the algorithm.

• Note that the CESampler abstract class implements the MVSampleIfc and RNStreamControlIfc interfaces. Thus, it can act as a multivariate sampler and can control the underlying random number stream. The four most relevant functions are:

• updateParameters(elites: List<DoubleArray>) - This function should update the parameters based on the points from the elite sample (held in a list of double arrays).

• parameters(): DoubleArray - This function returns the current parameter values.

• initializeParameters(values: DoubleArray) - This function is used to specify the initial values of the parameters.

• hasConverged(): Boolean - This function should return true if the parameters of the underlying sampling distribution are considered converged.

Initializing the Sampler

• For the CENormalSampler class, the implementation uses a randomly generated point in the solution space as the specification of the means of the multivariate normal distribution.
• Then, the input range of each point is used to approximate the initial values of the standard deviations for the distributions.
• The standard approximation for the standard deviation as a function of the range (\(R\)) is used: \(\sigma = R/4\).
• Thus, input variables that have a wider range of possible values start off with a higher variance parameter.
• The means associated with the normal distribution will converge to the recommended solution, while the standard deviations will continue to get closer and closer to zero as the process converges to a degenerate distribution.

Smoothing Constants and Convergence

• The CENormalSampler class requires the specification of the smoothing constant, \(\alpha\), for both the mean and the standard deviation parameters.
  • The default smoothing constants are set to \(0.85\).
• The user also has the ability to specify more or less variability via the initialVariabilityFactor property, which will be used to multiply the approximated initial values for the standard deviation parameters.
  • Its default value is 1.0.
• The user must also specify the cvThreshold parameter value.
  • This threshold is used to test if the parameters have converged.
  • The default implementation computes the coefficient of variation of the parameters and if the coefficient of variation goes below this threshold, then the parameter is consider to have converged.
  • If all the parameters have converged, then the hasConverged() function will return true.
  • The default coefficient of variation threshold is set to \(0.03\).

KSL R-SPLINE Main Iteration Loop

    override fun mainIteration() {
        // The initial (current) solution is randomly selected or specified by the user.
        // It will be the current solution until beaten by the SPLINE search process.
        // Call SPLINE for the next solution using the current sample size (m_k) and
        // current SPLINE oracle call limit (b_k).
        val splineSolution = spline(
            currentSolution,
            rSPLINESampleSize, splineCallLimit
        )
        currentSolution = splineSolution
        // capture the last solution
        solutionChecker.captureSolution(currentSolution)
    }

Overview of the SPLINE Function

• The SPLINE function is fairly complex.
• Given a specification of the sample size for the iteration, \(m_k\), the supplied solution may need additional replications evaluated, which will result in a new estimated solution.
• Then, an inner loop of the spline searching process commences.
• This involves two phases.
  • First, a search of a piece-wise linear interpolation of the objective function and then a search of the neighborhood resulting from the piece-wise linear interpolation search.
  • The details of this process can be found as algorithm 2 in Wang et al (2013)

Random Restart Solvers

• In general, optimization algorithms rely on a starting point to initialize the search process.
• The performance of the algorithm may be sensitive to the choice of the starting point, which may cause the search process to converge to a local optimal or to limit the exploration of the search space.
• A random restart solver is a solver that starts another solver using different randomly generated starting points.

class RandomRestartSolver(
    val restartingSolver: Solver,
    maxNumRestarts: Int = defaultMaxRestarts,
    streamNum: Int = 0,
    streamProvider: RNStreamProviderIfc = KSLRandom.DefaultRNStreamProvider,
    name: String? = null
) : StochasticSolver(
    restartingSolver.problemDefinition, restartingSolver.evaluator, maxNumRestarts,
    restartingSolver.replicationsPerEvaluation, streamNum, streamProvider, name
)

Random Restart Solver Main Iteration

• In sequential mode (the default), one inner solver is reused, run to completion once per restart.
• Since we want independent solutions, we clear a solution cache if it is being used.
• A new starting point is generated and assigned; the solver runs all its iterations (or until stopping). The best solution is captured.

    private fun sequentialMainIteration() {
        // clear the evaluator cache between randomized runs, but allow caching during the run itself
        // this will cause new replications to be generated
        if (clearCacheBetweenRuns){
            evaluator.cache?.clear()
        }
        // randomly assign a new starting point
        val startPoint = startingPoint()
        restartingSolver.startingPoint = startPoint
        // run the solver until it finds a solution
        restartingSolver.runAllIterations()
        numOracleCalls = numOracleCalls + restartingSolver.numOracleCalls
        numReplicationsRequested = numReplicationsRequested + restartingSolver.numReplicationsRequested
        // get the best solution from the solver run
        val bestSolution = restartingSolver.bestSolution
        currentSolution = bestSolution
    }

Applying a Solver

• To use a solver, we need at a minimum the following:
  • an instance of the ProblemDefinition class
  • an instance of something that implements the ModelBuilderIfc interface
  • an instance of the EvaluatorIfc interface
• The companion object of the Solver class provides factory functions that facilitate the creation of the available solvers with default configurations.
  • The factory functions will create an evaluator for executing the model that is configured to run the simulations locally (within the same execution thread) and is configured to use a MemorySolutionCache() to store and reuse solutions requested by the solver.

Supplying a Model Builder

object BuildRQModel : ModelBuilderIfc {
    override fun build(
        modelConfiguration: Map<String, String>?,
        experimentRunParameters: ExperimentRunParametersIfc?
    ): Model {
        val reorderQty: Int = 2
        val reorderPoint: Int = 1
        val model = Model("RQInventoryModel")
        val rqModel = RQInventorySystem(model, reorderPoint, reorderQty, "Inventory")
        rqModel.initialOnHand = 0
        rqModel.demandGenerator.initialTimeBtwEvents = ExponentialRV(1.0 / 3.6)
        rqModel.leadTime.initialRandomSource = ConstantRV(0.5)
        model.lengthOfReplication = 20000.0
        model.lengthOfReplicationWarmUp = 10000.0
        model.numberOfReplications = 40
        return model
    }
}

Summary of Classe

• ProblemDefinition - This class facilitates the definition of an optimization problem involving an objective function, deterministic linear constraints, deterministic functional constraints, and response constraints.
• PenaltyFunction - This abstract class defines the interface for penalty functions that can be used to handle constraint violations within an optimization problem, bound to the specific constraint each instance serves.
• InputDefinition - This class facilitates the definition of input variables for an optimization problem.
• ConstraintIfc - This interface defines the behavior of deterministic constraints for an optimization problem. There are two types of deterministic constraints specified by the LinearConstraint and FunctionalConstraint classes.
• ResponseConstraint - This class models constraints in a simulation optimization problem that are based on responses from the execution of the simulation model.

Summary of Classe

• Solution - This fundamental class encapsulates the inputs, outputs, and constraints associated with an evaluation request placed by a solver.
• Solver - This is an abstract base class the represents the concepts needed by all simulation optimization algorithms.
• StochasticSolver - This is an abstract class that represents solvers that use randomness in the search process. It has concrete implementations of StochasticHillClimber, SimulatedAnnealing, CrossEntropySolver, RSplineSovler, and RandomRestartSolver.
• ReplicationPerEvaluationIfc - This interface facilitates algorithms that change the number of replications during the search process. It has concrete implementations of FixedReplicationsPerEvaluation and FixedGrowthRateReplicationSchedule.
• SolutionEqualtyIfc - This interface as its variations allow for different processes to compare or check for solution equality, which is essential in implementing stopping criteria for algorithms.