1.Calibration framework and implementation

Any simulator, SimMobility included, has a set of parameters. Fixing an input, the output of the simulator varies depending of the values of such parameters. These values can be considered "good" if, providing the same input both to the real system being simulated and the simulator, the output given by the simulator is 'close' to the output measured in the real system. Calibration is the process that is meant to provide these "good" values[1].

Our calibration algorithm is based on the Simultaneous Perturbation method [2]. The figure below outlines the overall calibration framework for SimMobility Short-Term (ST). It is composed of three steps: Initialization, Perturbation, and Update parameter. Briefly, we first generate an initial vector of parameter values. Then, we will update from step to step this vector. In particular, given the current vector at a certain iteration, we first generate a random vector, of the same length of the parameter vector, which we call perturbation. We first add it to the current vector, obtaining the plus perturbed vector, we set the parameters in the simulator accordingly and we run the simulation. Then, we repeat the same operation for the minus perturbed vector, which is obtained by subtracting the perturbation from the current vector (instead of adding). Based on the outputs obtained with the plus and minus perturbed vectors, we compute the parameter vector for the next iteration.

Detailed framework of SimMobility ShortTerm calibration (wspsaNTE_enhanced.m)

1. Initialization

MT (SimMobility Mid-term)

• MT provides a Trip Chains (TC) table as a result of daily activity schedules. TC0 is the initial input trip chain table that has to be calibrated. Then, we should convert this trip chain table into Origin-Destination (OD) vector format.

Demand mapping (PostF_inLoop.R)

• Role:
  • F-function is a demand mapping algorithm that converts the TC into OD.
  • Also, F-function truncates the demand based on the capacity constraint.
  • Note that, Ffunction.R produces ODdata.csv and theta_OD.csv. This Ffunction.R is required to be implemented only in the case of regeneration of Trip Chain table by MT.
• Process: ODdata.csv, theta_OD.csv from Ffunction.R is the input for PostF_inLoop.R. With the predetermined Simulation start/end time and capacity constraints, we sum all trips in the same time-dependent origin/destination to produce Truncated_theta_od.csv.
• Note that:
  • PostF_inLoop.R: For initial run
  • PostF_inLoop_Plus.R: For the plus perturbed vector
  • PostF_inLoop_Minus.R: For the minus perturbed vector
  • PostF_inLoop_SPSA.R: For I of W-matrix

Configuration file (Configure_Whole_set.m)

• Role:
  • Configuration file includes the following configurations:
    • Driving behavior (OD) parameters: Initial value, Upper/Lower bound
    • Route choice (RC) parameters: Initial value, Upper/Lower bound
    • OD parameters (OD): Upper/Lower bound only
    • Algorithm parameters:
      • Alpha = 0.602
      • Convergence condition
      • Others: Free-flow interval / Warm up time / Number of replication
• Process: Read by the wspsaNTE_enhanced.m once initiated. Then, it provides the vector for initial parameter (θ_0): ‘Initial param’

where N_OD, N_DB and N_RC are the number of origin-destination parameters, driving behavior parameters and route choice parameters, respectively.

Running simulation (CHANGEPARAMETERS_pf.m)

• At each iteration, in order to find the parameter vector of the next iteration, we run SimMobiltiy Short Term three times. The CHANGEPARAMETERS_pf.m is able to run these three simulations simulations simultaneously. The process is illustrated in Figure 2 and can be described as follows:

Component of CHANGEPARAMETERS_pf.m

Example of OD parameter change over one iteration

Example of trip chain change over one iteration

Evaluation function (FUNC.m)

• Objective value
  • FUNC.m calculates the ‘Objective value’ for each measurements including counts, travel-time, and a-priori values.
  • We define the GLS objective function by giving weight with variability over different observation days of measurement. Then, the hybrid objective function is given as:

• Objective vector
  • Objective vector for each measurements are:

* These vectors are used in the gradient approximation. 

• Goodness of fit
  • RMSN measures the performance of the estimator with the average magnitude of the error. Because of the pre-squared error term, the index yields relatively high weights to large errors. This RMSN seems to be appropriate to deal with the large errors in calibration problem (compared with the others (e.g., MAE) which measure the accuracy of estimators with the individual differences with the equal weight in the average).
  • For count, travel-time, and a-priori, we can define the RMSN as:

2. Perturbation

Normalization of parameters

• Parameter vector
  • We define parameter vector for k-th iteration as:

• Normalized parameter vector
  • For each component i, we normalize with upper / lower boundary as:

Perturbations

• The random perturbation vector ∆_k is generated. Each component ∆_ki can be either {-1,1} with equal probability. Each component is extracted independently.
• Then, the plus and minus normalized perturbed vectors are:

* To remove normalization, these vectors are transformed into:

Simulation run with perturbed parameters

• After checking whether the two perturbed parameter vectors respect the constraints, we run Short Term (ST) twice (one with the plus perturbed vector and one with the minus perturbed vector) though CHANGEPARAMETERS_pf.m on parallel way using parpool as:

noPMrun=2;
    parpool('local', noPMrun);
    % Run SimMobility twice 
    parfor pf=1:noPMrun % pf=1: Plus, pf=2: Minus
        if pf==1 % Plus perturbed vector 
            dbparams =  sim_plus((end-no_rcparams-no_dbparam+1):(end-no_rcparams));
            rcparams =  sim_plus((end-no_rcparams+1):end);
        end
        if pf==2 % Minus perturbed vector
            dbparams =  sim_minus((end-no_rcparams-no_dbparam+1):(end-no_rcparams));
            rcparams =  sim_minus((end-no_rcparams+1):end);
        end
        PMstructure(pf) = CHANGEPARAMETERS_pf(dbparams, rcparams, iterations, pf);
    end
    delete(gcp);

Evaluate simulation result with perturbed parameters from plus / minus perturbed vector

• We evaluate FUNC.m to get the objective vector from plus and minus perturbed vector respectively.
• Note that CHANGEPARAMETERS_pf.m and FUNC.m have been thoroughly described in previous section.

3. Update parameters

The parameter vector update, to produce the parameter vector of the next iteration, is based on a gradient approximation vector, which is computed as follows.

Gradient approximation calculation

• In order to compute the gradient approximation, we first introduce the W matrix:

and define W_i its _i_-th row.

• The i-th component of the gradient approximation vector at iteration k is:

• The numerator is a vector that represents the gap between the two objective vectors, obtained with plus and minus perturbed parameter vectors, and can be expressed as:

• Finally, the gradient will be given as:

Update parameters

The value of the gradient approximation above is one of the possible vectors that we can use in SimMobility Short Term to update the parameter vector.

• It is possible to use the following types of vectors in the update rule:
  • g ̂_k: Weighted gradient as derived above (pf==1)
  • -g ̂_k: Opposite direction of the weighted gradient g ̂_k (pf==2)
  • And 〖g ̂_k〗_spsa: Weight matrix defined as I matrix for normal SPSA (pf==3)
• Then, the parameter vector is updated as:

Simulation run with updated parameters

• After checking the parameter constraints, we run ST twice though CHANGEPARAMETERS_pf.m on parallel way using parpool as:

noPMrun=3;
     parpool('local', noPMrun);
     parfor pf =1:noPMrun
      if pf == 1 % Plus perturbed vector (with Original gradient)
          dbparams_sim =  simulation((end-no_rcparams-no_dbparam+1):(end-no_rcparams)); 
          rcparams_sim =  simulation((end-no_rcparams+1):end);       
      end
      if pf == 2 % Minus perturbed vector (with Inverse gradient)
         dbparams_sim =  simulation_Inverse((end-no_rcparams-no_dbparam+1):(end-no_rcparams)); 
         rcparams_sim =  simulation_Inverse((end-no_rcparams+1):end);       
      end     
      if pf == 3 % SPSA side
         dbparams_sim =  simulation_0((end-no_rcparams-no_dbparam+1):(end-no_rcparams)); 
         rcparams_sim =  simulation_0((end-no_rcparams+1):end);       
      end   
         PMstructure(pf) = CHANGEPARAMETERS_pf(dbparams_sim, rcparams_sim, iterations, pf);
     end 
     delete(gcp);

Evaluate simulation result with three updated parameter

• We evaluate FUNC.m to get the objective value from the three updated parameter vectors above, Z_({pf==1}), Z_({pf==2}),and Z_({pf==3}).

Select the best parameter vector

• We take the best parameter vector among following 5 parameter vectors (θ(+), θ(-), θ ̂_(k+1_({pf==1}) ),θ ̂_(k+1_({pf==2}) ),θ ̂_(k+1_({pf==3}) )), which yields:

###References:

• [1] AA.VV. (n.d.). Fundamentals of Traffic Simulation. (J. Barcelo, Ed.). Springer.
• [2] Spall, J.C., 1998. An overview of the simultaneous perturbation method for efficient optimization. Johns Hopkins APL Tec. Dig. 19, 482–492.