在 R 中为自定义函数运行遗传算法时出错 IDM 的功能使用一组参数运行函数

问题描述

目标

我想估计智能驾驶员跟车模型 (IDM) 的“最佳”参数。 “最佳”是指那些在观察速度和预测速度之间产生最小均方根误差的参数。下面显示一个可重复的示例,我成功地使用网格搜索找到了最佳参数,但在运行遗传算法时却没有成功。

IDM 的功能

R 中的以下 IDM 函数接受 6 个参数并输出 3 列的数据帧,加速度 a_i,速度 v_i 和距离 g_x_i

calculate_IDM <- function(A_i,v_0,g_0,g_t_i,b_i,small_delta){
  
  ####################
  ## Allocate Vectors
  ####################
  
  # acceleration rate
  a_i <- rep(NA_real_,time_length)
  
  # speed
  v_i <- rep(NA_real_,time_length)
  
  
  # position
  x_i <- rep(NA_real_,time_length)
  
  # spacing
  g_x_i <- rep(NA_real_,time_length)
  
  # speed difference
  delta_v_i <- rep(NA_real_,time_length)
  
  
  # desired spacing
  g_star_i <- rep(NA_real_,time_length)
  
  
  
   ##########################################
  ## Initial values for Following vehicle
  ##########################################
  
  # speed
  v_i[1] <- v_i_first
  
  # position
  x_i[1] <- x_i_first
  
  # spacing
  g_x_i[1] <- xn1_first - x_i_first
  
  # speed difference
  delta_v_i[1] <- v_i_first - vn1_first
  
  # desired spacing
  g_star_i[1] <- g_0 + max(0,(v_i[1] * g_t_i) + ((v_i[1] * delta_v_i[1]) / (2 * sqrt((A_i * b_i)))))
  
  
  # acceleration rate
  a_i[1] <- A_i * (1 - ((v_i[1] / v_0)^small_delta) - ((g_star_i[1] / g_x_i[1])^2))
  
  # a_i[1] <- ifelse(is.nan(a_i[1]),A_i,a_i[1])
  
  
  
  # speed
  v_i[2] <- v_i[1] + (a_i[1] * time_frame)
  
  ### if the speed is negative,make it zero
  v_i[2] <- ifelse(v_i[2] < 0,v_i[2])
  
  
  
  # position
  x_i[2] <- x_i[1] + (v_i[1] * time_frame) + (0.5 * a_i[1] * (time_frame)^2)
  
  # spacing
  g_x_i[2] <- xn1_complete[2] - x_i[2]
  
  
  # speed difference
  delta_v_i[2] <- v_i[2] - vn1_complete[2]
  
  
  
  ####################
  ## IDM Calculations
  ####################
  
  
  for (t in 2:(time_length-1)) { 
    
    # desired spacing
    g_star_i[t] <- g_0 + max(0,(v_i[t] * g_t_i) + ((v_i[t] * delta_v_i[t]) / (2 * sqrt((A_i * b_i)))))
    
    
    # acceleration rate
    a_i[t] <- A_i * (1 - ((v_i[t] / v_0)^small_delta) - ((g_star_i[t] / g_x_i[t])^2))
    
    # a_i[t] <- ifelse(is.nan(a_i[t]),a_i[t])
    
    
    
    # speed
    v_i[t+1] <- v_i[t] + (a_i[t] * time_frame)
    
    ### if the speed is negative,make it zero
    v_i[t+1] <- ifelse(v_i[t+1] < 0,v_i[t+1])
    
    
    
    # position
    x_i[t+1] <- x_i[t] + (v_i[t] * time_frame) + (0.5 * a_i[t] * (time_frame)^2)
    
    # spacing
    g_x_i[t+1] <- xn1_complete[t+1] - x_i[t+1]
    
    
    # speed difference
    delta_v_i[t+1] <- v_i[t+1] - vn1_complete[t+1]
    
    
  }
  
  data.frame(a_i,v_i,g_x_i)
}

使用一组参数运行函数

要运行上述函数,需要先导车辆的速度和时间向量:

# Time
last_time <- 300 ## s
time_frame <- 0.1 ## s
Time <- seq(from = 0,to = last_time,by = time_frame)
time_length <- length(Time)


v_i_first <- 0
x_i_first <- 5


## Lead vehicle
vn1_first <- 0 ## first speed m/s
xn1_first <- 100 ## position of lead vehicle front center m
bn1_complete <- c(rep(x = 4,length.out = time_length-2951),rep(x = 0,length.out = time_length-50)) ## acceleration rate 



#############################################
### Complete speed trajectory of Lead vehicle
#############################################

vn1_complete <- rep(NA_real_,time_length) ### an empty vector 
xn1_complete <- rep(NA_real_,time_length) ### an empty vector 

vn1_complete[1] <- vn1_first
xn1_complete[1] <- xn1_first

for (t in 2:time_length) { 
  
  ### Lead vehicle calculations
  vn1_complete[t] <- vn1_complete[t-1] + (bn1_complete[t-1] * time_frame)
  
  
  xn1_complete[t] <- xn1_complete[t-1] + (vn1_complete[t-1] * time_frame) + (0.5 * bn1_complete[t-1] * (time_frame)^2)
  
}

现在,我可以应用这个函数了:

## one example
dff <- calculate_IDM(4,30,6.5,1,4,2)
head(dff)
       a_i       v_i    g_x_i
1 3.981274 0.0000000 95.00000
2 3.978206 0.3981274 95.00009
3 3.973594 0.7959480 95.00039
4 3.967446 1.1933075 95.00093
5 3.959771 1.5900521 95.00176
6 3.950581 1.9860292 95.00296

使用网格搜索找到最佳参数:

观察到的速度和参数列表如下:

library(tidyverse)

obs_data <- tibble(
obs_v_i = dff$v_i
)

# Parameters list
parameters_grid <- list(
  A_i = c(2,4),v_0 = c(27,30),g_0 = c(6.5,7),g_t_i = c(1,3),b_i = c(4,5),small_delta = c(2,3)
) %>% 
  expand.grid() %>% 
  as_tibble()

适应度函数和两个例子如下:

# fitness function
fitness_func <- function(obs_data,small_delta) {
  
  df <- cbind(obs_data,calculate_IDM(A_i,small_delta))
  
  sqrt(sum((df$obs_v_i - df$v_i)^2)/nrow(df))
  
}

> fitness_func(obs_data,2)
[1] 0
> fitness_func(obs_data,2,27,7,3,5,3)
[1] 1.406937

现在我可以使用 rowwise() 中的 dplyr 函数进行网格搜索

parameters_grid %>% 
  rowwise() %>% 
  mutate(RMSE = fitness_func(obs_data,small_delta))

# A tibble: 64 x 7
# Rowwise: 
     A_i   v_0   g_0 g_t_i   b_i small_delta    RMSE
   <dbl> <dbl> <dbl> <dbl> <dbl>       <dbl>   <dbl>
 1     2    27   6.5     1     4           2 1.68   
 2     4    27   6.5     1     4           2 0.213  
 3     2    30   6.5     1     4           2 1.65   
 4     4    30   6.5     1     4           2 0      
 5     2    27   7       1     4           2 1.68   
 6     4    27   7       1     4           2 0.218  
 7     2    30   7       1     4           2 1.65   
 8     4    30   7       1     4           2 0.00794
 9     2    27   6.5     3     4           2 1.57   
10     4    27   6.5     3     4           2 0.814  
# ... with 54 more rows

遗传算法错误

你可以想象,如果参数列表更大,它会显着增加计算时间。所以,我想运行遗传算法。使用示例 here,我尝试使用 GA 库来估计参数但出现错误

library(GA)
GA <- ga(type = "real-valued",fitness =  -fitness_func(obs_data,small_delta),lower = c(2,2),upper = c(4,popSize = 5,maxiter = 10,run = 10)

 Error in calculate_IDM(A_i,small_delta) : 
  object 'g_0' not found 

请让我知道我在这里做错了什么。

解决方法

fitness 中记录的 ?ga

适应度函数,任何允许的 R 函数,它将表示潜在解决方案的单个字符串作为输入,并返回一个描述其“适应度”的数值。

所以,我们可以将它包装成一个带有两个参数的函数,然后使用 fitness_func 参数作为 x[1],x[2],...,x[6] 这将是与 lowerupper 绑定值的长度相同。在这里,我们也可以单独传递data

library(GA)
GA <- ga(type = "real-valued",fitness =  function(dat,x) {-fitness_func(dat,x[1],x[2],x[3],x[4],x[5],x[6])},dat = obs_data,lower = c(2,27,6.5,1,4,2),upper = c(4,30,7,3,5,3),popSize = 5,maxiter = 1000,run = 100)
#GA | iter = 1 | Mean = -0.5668704 | Best = -0.3523867
#GA | iter = 2 | Mean = -0.3762976 | Best = -0.3523867
#GA | iter = 3 | Mean = -0.3529940 | Best = -0.3523867
#GA | iter = 4 | Mean = -0.3523867 | Best = -0.3523867
#GA | iter = 5 | Mean = -0.3523867 | Best = -0.3523867
#GA | iter = 6 | Mean = -0.3523867 | Best = -0.3523867
#GA | iter = 7 | Mean = -0.3523867 | Best = -0.3523867
#GA | iter = 8 | Mean = -0.3640060 | Best = -0.3523867
#...
#GA | iter = 519 | Mean = -0.08506463 | Best = -0.08505393
#GA | iter = 520 | Mean = -0.08505440 | Best = -0.08505393
#GA | iter = 521 | Mean = -0.14507196 | Best = -0.08505393
#GA | iter = 522 | Mean = -0.08505393 | Best = -0.08505393
#GA | iter = 523 | Mean = -0.08505393 | Best = -0.08505393
#GA | iter = 524 | Mean = -0.11238973 | Best = -0.08505393
#GA | iter = 525 | Mean = -0.31888465 | Best = -0.08505393
#GA | iter = 526 | Mean = -0.09641056 | Best = -0.08505393

虽然最后有警告