代码之家  ›  专栏  ›  技术社区  ›  umair durrani

在R中为自定义函数运行遗传算法时出错

  •  0
  • umair durrani  · 技术社区  · 4 年前

    目标

    我想估计智能驾驶员跟车模型(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, 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, 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, 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()
    

    健身功能和2个示例如下:

    # Fitness function
    fitness_func <- function(obs_data,
                             A_i, 
                             v_0, 
                             g_0,
                             g_t_i, 
                             b_i, 
                             small_delta) {
      
      df <- cbind(obs_data, calculate_IDM(A_i, 
                                    v_0, 
                                    g_0,
                                    g_t_i, 
                                    b_i, 
                                    small_delta))
      
      sqrt(sum((df$obs_v_i - df$v_i)^2)/nrow(df))
      
    }
    
    > fitness_func(obs_data, 4, 30, 6.5, 1, 4, 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,
                                 A_i, 
                                 v_0, 
                                 g_0,
                                 g_t_i, 
                                 b_i, 
                                 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,
                                     A_i, 
                                     v_0, 
                                     g_0,
                                     g_t_i, 
                                     b_i, 
                                     small_delta),
             lower = c(2, 27, 6.5, 1, 4, 2), upper = c(4, 30, 7, 3, 5, 3), 
             popSize = 5, maxiter = 10, run = 10)
    
     Error in calculate_IDM(A_i, v_0, g_0, g_t_i, b_i, small_delta) : 
      object 'g_0' not found 
    

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

    0 回复  |  直到 4 年前
        1
  •  2
  •   akrun    4 年前

    这个 fitness 如文件所述 ?ga

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

    所以,我们可以把它包装成一个有两个参数的函数,然后使用 fitness_func 论点如 x[1] , x[2] , ..., x[6] 其长度与 lower upper 约束值。在这里,我们也可以通过 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
    

    但最后有一个警告