Simulated Annealing: Global Optimization Without Gradients

Part 36 of 12-Week BusinessMath Series

What You’ll Learn

• Understanding simulated annealing for global optimization
• Cooling schedules: exponential, linear, adaptive
• Acceptance probability and the Metropolis criterion
• When to use simulated annealing vs. gradient methods
• Escaping local minima through controlled randomness
• Parameter tuning: initial temperature, cooling rate, iterations

The Problem

Many business optimization problems have multiple local minima:

• Production scheduling with setup costs (discontinuous objective)
• Portfolio optimization with transaction costs and lot sizes
• Facility location with discrete choices (city A vs. city B)
• Hyperparameter tuning for machine learning models

Gradient-based methods get stuck in local minima. Need global search capability.

The Solution

Simulated annealing mimics the physical process of metal cooling: accept worse solutions with probability that decreases over time. This allows escaping local minima early while converging to global optimum later.

Pattern 1: Portfolio Optimization with Transaction Costs

Business Problem: Rebalance portfolio from suboptimal starting point, but minimize transaction costs (non-smooth objective).

import BusinessMath

// Portfolio with transaction costs (20 assets for clarity)
let numAssets = 20

// Create a suboptimal starting portfolio: heavily concentrated in first 5 assets
// These happen to be LOW return assets - clear opportunity to improve!
var currentWeights = [Double](repeating: 0.0, count: numAssets)
// First 5 assets: 60% of portfolio (12% each) - these are low-return!
for i in 0..<5 {
    currentWeights[i] = 0.12
}
// Remaining 15 assets: 40% of portfolio (2.67% each) - these include high-return!
for i in 5..
            
              = 5 && i < 15 && j >= 5 && j < 15) || (i >= 15 && j >= 15) return sameTier ? Double.random(in: 0.5...0.7) : Double.random(in: 0.2...0.4) } } func portfolioObjective(_ weights: VectorN
              
                ) -> Double { // 1. Portfolio expected return (we want to MAXIMIZE this) var portfolioReturn = 0.0 for i in 0..
                
                  >( config: config, searchSpace: searchSpace ) // Create initial guess from current weights let initialGuess = VectorN(currentWeights) // Define constraints (FIXED API - no 'constraint:' or 'tolerance:' parameters!) let constraints: [MultivariateConstraint
                  
                    >] = [ // Equality: Sum to 1 (must be fully invested) .equality { weights in (0..
                    
                       abs($1.change) }.prefix(8) print("\nTop 8 Position Changes:") print("───────────────────────────────────────────────────────────") print("Asset Return Old Weight New Weight Change Action") print("───────────────────────────────────────────────────────────") for change in changes { let direction = change.change > 0 ? "BUY " : "SELL" print(String(format: " %2d %5.2f%% %6.2f%% %6.2f%% %+6.2f%% %s", change.index, change.return * 100, change.oldWeight * 100, change.newWeight * 100, change.change * 100, direction)) } print("═══════════════════════════════════════════════════════════") print("\n💡 Key Insight:") print(" The optimizer balanced improving returns (moving to high-return assets)") print(" with minimizing transaction costs. Notice it didn't eliminate all") print(" low-return holdings - the transaction costs made that too expensive.") 
                    
                  
                
              
            

Pattern 2: Configuration Comparison

Pattern: Compare different cooling configurations.

// Simple test function: Rastrigin in 2D (many local minima)
func rastrigin(_ x: VectorN
            
              ) -> Double { let A = 10.0 let n = 2.0 return A * n + (0..<2).reduce(0.0) { sum, i in sum + (x[i] * x[i] - A * cos(2 * .pi * x[i])) } } let searchSpace = [(-5.12, 5.12), (-5.12, 5.12)] // Test different configurations let configs: [(name: String, config: SimulatedAnnealingConfig)] = [ ("Fast", .fast), ("Default", .default), ("Thorough", .thorough), ("Custom (slow)", SimulatedAnnealingConfig( initialTemperature: 100.0, finalTemperature: 0.001, coolingRate: 0.98, // Slower cooling maxIterations: 20_000, perturbationScale: 0.1, reheatInterval: nil, reheatTemperature: nil, seed: 42 )) ] print("Configuration Comparison (Rastrigin Function)") print("═══════════════════════════════════════════════════════════") print("Config | Final Value | Iterations | Acceptance Rate") print("───────────────────────────────────────────────────────────") for (name, config) in configs { let optimizer = SimulatedAnnealing
              
                >( config: config, searchSpace: searchSpace ) let result = optimizer.optimizeDetailed( objective: rastrigin, initialSolution: VectorN([2.0, 3.0]) ) let rate = result.acceptanceRate * 100 print("\(name.padding(toLength: 15, withPad: " ", startingAt: 0)) | " + "\(result.fitness.formatted(.number.precision(.fractionLength(6))).padding(toLength: 11, withPad: " ", startingAt: 0)) | " + "\(String(format: "%10d", result.iterations)) | " + "\(rate.formatted(.number.precision(.fractionLength(1))))%") } print("\nRecommendation: Use .default for most problems, .thorough for difficult landscapes") 
              
            

Pattern 3: Ackley Function (Multimodal Optimization)

Pattern: Global optimization for highly multimodal functions.

// Ackley function: highly multimodal with many local minima
// Global minimum at (0, 0) with value 0
func ackley(_ x: VectorN
            
              ) -> Double { let a = 20.0 let b = 0.2 let c = 2.0 * .pi let d = 2 // dimensions let sum1 = (0..
              
                >( config: .thorough, searchSpace: searchSpace ) print("Ackley Function Optimization (Highly Multimodal)") print("═══════════════════════════════════════════════════════════") // Start from a poor initial guess let initialGuess = VectorN([4.0, -3.5]) let result = sa.optimizeDetailed( objective: ackley, initialSolution: initialGuess ) print("Optimization Results:") print(" Solution: (\(result.solution[0].formatted(.number.precision(.fractionLength(4)))), " + "\(result.solution[1].formatted(.number.precision(.fractionLength(4)))))") print(" Function Value: \(result.fitness.formatted(.number.precision(.fractionLength(6)))) (target: 0.0)") print(" Iterations: \(result.iterations)") print(" Acceptance Rate: \((result.acceptanceRate * 100).formatted(.number.precision(.fractionLength(1))))%") print(" Converged: \(result.converged)") print(" Reason: \(result.convergenceReason)") // Distance from global optimum let distanceFromOptimum = sqrt(result.solution[0]*result.solution[0] + result.solution[1]*result.solution[1]) print(" Distance from global optimum: \(distanceFromOptimum.formatted(.number.precision(.fractionLength(4))))") 
              
            

How It Works

Simulated Annealing Algorithm

1. Initialize: Set T = T_0, x = x_0
2. Generate Neighbor: x’ = random perturbation of x
3. Calculate ΔE: ΔE = f(x’) - f(x)
4. Accept/Reject:
   • If ΔE < 0: Always accept (improvement)
   • Else: Accept with probability P = exp(-ΔE / T)
5. Cool Down: T = α * T (exponential) or T = T - β (linear)
6. Repeat: Until T < T_min or max iterations

Acceptance Probability

Metropolis Criterion: P(accept worse solution) = exp(-ΔE / T)

Insight: Early (high T), accept almost anything. Late (low T), only accept improvements.

Cooling Schedule Impact

Problem: 100-variable portfolio optimization

Tradeoff: Slower cooling = better solution, more time

Real-World Application

Manufacturing: Batch Sizing with Setup Costs

Company: Electronics manufacturer optimizing production batch sizes for 15 products Challenge: Minimize total costs (inventory + setup) subject to demand and capacity

Problem Characteristics:

• Continuous variables: Batch sizes (treated as continuous for optimization)
• Setup costs: Fixed cost per production run (non-smooth)
• Holding costs: Penalty for excess inventory (smooth)
• Multiple local minima: ~100+ feasible configurations

Why Simulated Annealing:

• Setup costs create discontinuous objective
• Global search needed (many local minima)
• Can escape poor local solutions

Implementation:

import BusinessMath

import BusinessMath

// Model with 15 products
let numProducts = 15

// Weekly demand per product (units/week)
let demand: [Double] = [
    100.0, 150.0, 80.0, 200.0, 120.0,  // Products 0-4
    90.0, 175.0, 110.0, 140.0, 95.0,   // Products 5-9
    160.0, 85.0, 130.0, 105.0, 145.0   // Products 10-14
]

// Fixed setup cost per production run ($)
let setupCost: [Double] = [
    500.0, 750.0, 400.0, 850.0, 600.0,  // Products 0-4
    450.0, 800.0, 550.0, 700.0, 480.0,  // Products 5-9
    720.0, 420.0, 650.0, 520.0, 680.0   // Products 10-14
]

// Holding cost per unit per week ($/unit/week)
let holdingCost: [Double] = [
    2.0, 3.5, 1.8, 4.0, 2.5,    // Products 0-4
    1.9, 3.8, 2.2, 3.0, 2.0,    // Products 5-9
    3.3, 1.7, 2.8, 2.1, 3.2     // Products 10-14
]

func totalCost(_ batchSizes: VectorN
            
              ) -> Double { var cost = 0.0 for i in 0..
              
                >( config: config, searchSpace: searchSpace ) print("\nRunning Simulated Annealing...") let result = sa.optimizeDetailed( objective: totalCost, initialSolution: naiveBatches ) print("\nOptimization Results:") print(" Converged: \(result.converged)") print(" Iterations: \(result.iterations)") print(" Final Temperature: \(result.finalTemperature.formatted(.number.precision(.fractionLength(4))))") print(" Acceptance Rate: \(result.acceptanceRate.percent(1))") let optimizedCost = result.fitness let costReduction = naiveCost - optimizedCost let percentReduction = (costReduction / naiveCost) * 100 let weeklySavings = costReduction print("\n═══════════════════════════════════════════════════════════") print("Cost Comparison:") print("═══════════════════════════════════════════════════════════") print("Naive Approach: $\(naiveCost.formatted(.number.precision(.fractionLength(2))))") print("Optimized (SA): $\(optimizedCost.formatted(.number.precision(.fractionLength(2))))") print("Cost Reduction: $\(costReduction.formatted(.number.precision(.fractionLength(2)))) (\(percentReduction.formatted(.number.precision(.fractionLength(1))))%)") print("Weekly Savings: $\(weeklySavings.formatted(.number.precision(.fractionLength(2))))") print("═══════════════════════════════════════════════════════════") // Show optimal batch sizes for top 5 products by demand let productInfo = (0..
                
                   $1.demand } print("\nOptimal Batch Sizes (Top 5 by Demand):") print("───────────────────────────────────────────────────────────") print("Product Demand Optimal Batch Runs/Week Setup $ Hold $") print("───────────────────────────────────────────────────────────") for info in productInfo.prefix(5) { let runsPerWeek = info.demand / info.optimalBatch let setupCostWeekly = runsPerWeek * info.setupCost let holdCostWeekly = (info.optimalBatch / 2.0) * info.holdingCost print(String(format: " %2d %5.0f %6.1f %5.2f $%5.0f $%5.0f", info.id, info.demand, info.optimalBatch, runsPerWeek, setupCostWeekly, holdCostWeekly)) } print("═══════════════════════════════════════════════════════════") 
                
              
            

Typical Results:

• Cost reduction: 10-15% vs. naive approach (batch size = demand)
• Weekly savings: $6,000-$10,000 depending on product mix
• Computation time: 30-60 seconds (acceptable for weekly planning)
• Solution quality: Near-optimal (escapes local minima that gradient methods miss)

Try It Yourself

import Foundation
import BusinessMath

// Portfolio with transaction costs (20 assets for clarity)
let numAssets = 20

// Create a suboptimal starting portfolio: heavily concentrated in first 5 assets
// These happen to be LOW return assets - clear opportunity to improve!
var currentWeights = [Double](repeating: 0.0, count: numAssets)
// First 5 assets: 60% of portfolio (12% each) - these are low-return!
for i in 0..<5 {
    currentWeights[i] = 0.12
}
// Remaining 15 assets: 40% of portfolio (2.67% each) - these include high-return!
for i in 5..
              
                = 5 && i < 15 && j >= 5 && j < 15) || (i >= 15 && j >= 15) return sameTier ? Double.random(in: 0.5...0.7) : Double.random(in: 0.2...0.4) } } @MainActor func portfolioObjective(_ weights: VectorN
                
                  ) -> Double { // 1. Portfolio expected return (we want to MAXIMIZE this) var portfolioReturn = 0.0 for i in 0..
                  
                    >( config: config, searchSpace: searchSpace ) // Create initial guess from current weights let initialGuess = VectorN(currentWeights) // Define constraints let constraints: [MultivariateConstraint
                    
                      >] = [ // Equality: Sum to 1 (must be fully invested) .equality { weights in (0..
                      
                         abs($1.change) }.prefix(8) print("\nTop 8 Position Changes:") print("───────────────────────────────────────────────────────────") print("Asset Return Old Weight New Weight Change Action") print("───────────────────────────────────────────────────────────") for change in changes { let direction = change.change > 0 ? "BUY " : "SELL" print("\("\(change.index)".paddingLeft(toLength: 4))\(change.return.percent().paddingLeft(toLength: 9))\(change.oldWeight.percent(2).paddingLeft(toLength: 13))\(change.newWeight.percent().paddingLeft(toLength: 12))\(change.change.percent(2).paddingLeft(toLength: 8))\(direction.paddingLeft(toLength: 12))") } print("═══════════════════════════════════════════════════════════") print("\n💡 Key Insight:") print(" The optimizer balanced improving returns (moving to high-return assets)") print(" with minimizing transaction costs. Notice it didn't eliminate all") print(" low-return holdings - the transaction costs made that too expensive.") // MARK: - Configuration Comparison // Simple test function: Rastrigin in 2D (many local minima) func rastrigin(_ x: VectorN
                        
                          ) -> Double { let A = 10.0 let n = 2.0 return A * n + (0..<2).reduce(0.0) { sum, i in sum + (x[i] * x[i] - A * cos(2 * .pi * x[i])) } } let searchSpace_config = [(-5.12, 5.12), (-5.12, 5.12)] // Test different configurations let configs: [(name: String, config: SimulatedAnnealingConfig)] = [ ("Fast", .fast), ("Default", .default), ("Thorough", .thorough), ("Custom (slow)", SimulatedAnnealingConfig( initialTemperature: 100.0, finalTemperature: 0.001, coolingRate: 0.98, // Slower cooling maxIterations: 20_000, perturbationScale: 0.1, reheatInterval: nil, reheatTemperature: nil, seed: 42 )) ] print("Configuration Comparison (Rastrigin Function)") print("═══════════════════════════════════════════════════════════") print("Config | Final Value | Iterations | Acceptance Rate") print("───────────────────────────────────────────────────────────") for (name, config) in configs { let optimizer = SimulatedAnnealing
                          
                            >( config: config, searchSpace: searchSpace_config ) let result = optimizer.optimizeDetailed( objective: rastrigin, initialSolution: VectorN([2.0, 3.0]) ) let rate = result.acceptanceRate print("\(name.padding(toLength: 15, withPad: " ", startingAt: 0)) | " + "\(result.fitness.formatted(.number.precision(.fractionLength(6))).padding(toLength: 11, withPad: " ", startingAt: 0)) | " + "\(String(format: "%10d", result.iterations)) | " + "\(rate.percent(1))") } print("\nRecommendation: Use .default for most problems, .thorough for difficult landscapes") // MARK: - Ackley Function (Multimodal Optimization) // Ackley function: highly multimodal with many local minima // Global minimum at (0, 0) with value 0 func ackley(_ x: VectorN
                            
                              ) -> Double { let a = 20.0 let b = 0.2 let c = 2.0 * .pi let d = 2 // dimensions let sum1 = (0..
                              
                                >( config: .thorough, searchSpace: searchSpace_ackley ) print("Ackley Function Optimization (Highly Multimodal)") print("═══════════════════════════════════════════════════════════") // Start from a poor initial guess let initialGuess_ackley = VectorN([4.0, -3.5]) let result_ackley = sa_ackley.optimizeDetailed( objective: ackley, initialSolution: initialGuess_ackley ) print("Optimization Results:") print(" Solution: (\(result_ackley.solution[0].number(4)), " + "\(result_ackley.solution[1].number(4)))") print(" Function Value: \(result_ackley.fitness.number(6)) (target: 0.0)") print(" Iterations: \(result_ackley.iterations)") print(" Acceptance Rate: \(result_ackley.acceptanceRate.percent(1))") print(" Converged: \(result_ackley.converged)") print(" Reason: \(result_ackley.convergenceReason)") // Distance from global optimum let distanceFromOptimum = sqrt(result_ackley.solution[0]*result_ackley.solution[0] + result_ackley.solution[1]*result_ackley.solution[1]) print(" Distance from global optimum: \(distanceFromOptimum.number(4))") // MARK: - Real-World Example: Manufacturing Batch Sizing with Setup Costs print("\n\nManufacturing Batch Sizing Optimization") print("═══════════════════════════════════════════════════════════") print("\nProblem: 15 products with setup costs and inventory holding costs") // Model with 15 products let numProducts_batch = 15 // Weekly demand per product (units/week) let demand_batch: [Double] = [ 100.0, 150.0, 80.0, 200.0, 120.0, // Products 0-4 90.0, 175.0, 110.0, 140.0, 95.0, // Products 5-9 160.0, 85.0, 130.0, 105.0, 145.0 // Products 10-14 ] // Fixed setup cost per production run ($) let setupCost_batch: [Double] = [ 500.0, 750.0, 400.0, 850.0, 600.0, // Products 0-4 450.0, 800.0, 550.0, 700.0, 480.0, // Products 5-9 720.0, 420.0, 650.0, 520.0, 680.0 // Products 10-14 ] // Holding cost per unit per week ($/unit/week) let holdingCost_batch: [Double] = [ 2.0, 3.5, 1.8, 4.0, 2.5, // Products 0-4 1.9, 3.8, 2.2, 3.0, 2.0, // Products 5-9 3.3, 1.7, 2.8, 2.1, 3.2 // Products 10-14 ] func totalCost_batch(_ batchSizes: VectorN
                                
                                  ) -> Double { var cost = 0.0 for i in 0..
                                  
                                    >( config: config_batch, searchSpace: searchSpace_batch ) print("\nRunning Simulated Annealing...") let result_batch = sa_batch.optimizeDetailed( objective: totalCost_batch, initialSolution: naiveBatches ) print("\nOptimization Results:") print(" Converged: \(result_batch.converged)") print(" Iterations: \(result_batch.iterations)") print(" Final Temperature: \(result_batch.finalTemperature.number(4))") print(" Acceptance Rate: \(result_batch.acceptanceRate.percent(1))") let optimizedCost = result_batch.fitness let costReduction = naiveCost - optimizedCost let percentReduction = (costReduction / naiveCost) * 100 let weeklySavings = costReduction print("\n═══════════════════════════════════════════════════════════") print("Cost Comparison:") print("═══════════════════════════════════════════════════════════") print("Naive Approach: \(naiveCost.currency())") print("Optimized (SA): \(optimizedCost.currency())") print("Cost Reduction: \(costReduction.currency()) (\(percentReduction.number(1))%)") print("Weekly Savings: \(weeklySavings.currency())") print("═══════════════════════════════════════════════════════════") // Show optimal batch sizes for top 5 products by demand let productInfo = (0..
                                    
                                       $1.demand } print("\nOptimal Batch Sizes (Top 5 by Demand):") print("─────────────────────────────────────────────────────────────") print("Product Demand Optimal Batch Runs/Week Setup $ Hold $") print("─────────────────────────────────────────────────────────────") for info in productInfo.prefix(5) { let runsPerWeek = info.demand / info.optimalBatch let setupCostWeekly = runsPerWeek * info.setupCost let holdCostWeekly = (info.optimalBatch / 2.0) * info.holdingCost print("\("\(info.id)".paddingLeft(toLength: 7))\(info.demand.number(0).paddingLeft(toLength: 8))\(info.optimalBatch.number(1).paddingLeft(toLength: 15))\(runsPerWeek.number(2).paddingLeft(toLength: 11))\(setupCostWeekly.currency().paddingLeft(toLength: 9))\(holdCostWeekly.currency().paddingLeft(toLength: 11))") } print("═════════════════════════════════════════════════════════════") 
                                    
                                  
                                
                              
                            
                          
                        
                      
                    
                  
                
              

→ Full API Reference: BusinessMath Docs – Simulated Annealing Tutorial

Experiments to Try

1. Temperature Tuning: Test initial temperatures 0.1, 1.0, 10.0, 100.0
2. Cooling Rates: Compare α = 0.90, 0.95, 0.98, 0.99
3. Neighbor Generation: Different perturbation strategies for portfolio
4. Hybrid Approach: SA for global search, then local refinement with BFGS

Next Steps

Next Week: Week 11 begins with Nelder-Mead Simplex (another gradient-free method), then Particle Swarm Optimization, concluding with Case Study #5: Real-Time Portfolio Rebalancing.

Final Week: Week 12 covers reflections (What Worked, What Didn’t) and Case Study #6: Investment Strategy DSL.

Series: [Week 10 of 12] | Topic: [Part 5 - Advanced Methods] | Case Studies: [4/6 Complete]

Topics Covered: Simulated annealing • Global optimization • Cooling schedules • Metropolis criterion • Discrete optimization

Playgrounds: [Week 1-10 available] • [Next week: Nelder-Mead and particle swarm]