如何使用R的optim函数求解满足流量匹配要求的水面高程WSE

我需要确定水力参数水面高程（Water Surface Elevation, WSE），使得给定点位的河流计算流量与实测流量匹配。我有大量点位需要迭代处理，因此希望编写循环批量运行所有数据，目前先在指定的断面位置进行测试。
我期望optim()函数能找到最优WSE，将其代入其他水力公式（过流面积、湿周）估算河流流量。但目前optim()返回结果不符合预期，以下是测试代码和返回结果：

# 构造测试数据:
x <- seq(0, 50, length.out = 20)
z <- c(seq(100,60,length.out = 8), 58,56,56,58, seq(60,100, length.out = 8))
plot(x,z) # V型河道横断面

# 去除空值（真实数据中存在大量空值）:
x <- x[!is.na(x)]
z <- z[!is.na(z)]

# 曼宁系数:
N <- 0.05

# 目标流量（单位：立方英尺/秒）: 
Q <- 400

# 河道坡降（单位：ft/ft）
slope <- 0.004

# 构造函数最小化计算流量与实测流量的差值:
Q.calc <- function(WSE){
  
  # 计算过流面积:
  wet.area <- vector(length = length(x))
  for(i in 1:(length(x)-1)){
    if(max(z[i],z[i+1]) > WSE){ # 水面以上的点位赋值为0
      wet.area[i] <- 0
    } else {
      wet.area[i] <- (x[i+1] -x[i])*(WSE - (z[i+1] + z[i])/2)
    }
  }
  
  wet.area <- sum(wet.area) # 全断面过流面积求和
  
  # 计算湿周:
  wet.perim <- vector(length = length(z))
  for(i in 1:(length(z)-1)){
    if(max(z[i],z[i+1]) > WSE) {  # 水面以上的点位赋值为0
      wet.perim[i] <- 0
    } else {  # 计算水面以下部分的湿周
      wet.perim[i] <- sqrt(((x[i+1]-x[i])^2)+((abs(z[i+1]-z[i])^2)))
    }
  }
  wet.perim <- sum(wet.perim) # 全断面湿周求和
  
  # 计算水力半径:
  hyd.radius <- wet.area/wet.perim
  
  # 实际要优化的等式:
  # 计算Q.calc等于实测Q的条件:  
  (1.49/N)*wet.area*hyd.radius^(2/3)*slope^0.5 - Q # 优化时减去实测Q，让最优解下计算值等于实测值
}

optim(par = mean(z), 
      fn = Q.calc, 
      lower = min(z), 
      upper = max(z),
      method = 'Brent',
      control = list(maxit = 100000))

# 返回结果
$par
[1] 56

$value
[1] -400

$counts
function gradient 
      NA       NA 

$convergence
[1] 0

$message
NULL

• 当前定义的目标函数返回值为计算流量 - 目标流量Q，而optim()默认逻辑是寻找函数的最小值。当WSE等于断面最低点56时，过流面积为0，计算流量为0，返回值为0-400=-400，刚好是整个取值区间的最小值，所以optim()直接返回了WSE=56的结果，没有找到让计算流量等于目标流量的解。
• 你需要优化的目标应该是计算流量与目标流量差值的绝对值（或平方值），只有这类目标的最小值才会出现在差值为0的位置，符合你找WSE让计算流量匹配实测流量的需求。

只需要修改目标函数的返回值，改为返回差值的绝对值即可，同时增加湿周为0的异常处理，避免报错，修改后的完整代码如下：

# 构造测试数据
x <- seq(0, 50, length.out = 20)
z <- c(seq(100,60,length.out = 8), 58,56,56,58, seq(60,100, length.out = 8))
# 去除空值
x <- x[!is.na(x)]
z <- z[!is.na(z)]
# 曼宁系数
N <- 0.05
# 目标流量（立方英尺/秒）
Q <- 400
# 河道坡降
slope <- 0.004

# 修正后的目标函数，最小化计算流量与目标流量差值的绝对值
Q.calc <- function(WSE){
  # 计算过流面积
  wet.area <- vector(length = length(x))
  for(i in 1:(length(x)-1)){
    if(max(z[i],z[i+1]) > WSE){ 
      wet.area[i] <- 0
    } else {
      wet.area[i] <- (x[i+1] -x[i])*(WSE - (z[i+1] + z[i])/2)
    }
  }
  wet.area <- sum(wet.area) 
  # 计算湿周
  wet.perim <- vector(length = length(z))
  for(i in 1:(length(z)-1)){
    if(max(z[i],z[i+1]) > WSE) { 
      wet.perim[i] <- 0
    } else { 
      wet.perim[i] <- sqrt(((x[i+1]-x[i])^2)+((abs(z[i+1]-z[i])^2)))
    }
  }
  wet.perim <- sum(wet.perim) 
  # 避免湿周为0导致的除以0报错
  if(wet.perim == 0){
    hyd.radius <- 0
  } else {
    hyd.radius <- wet.area/wet.perim
  }
  # 计算流量，返回差值的绝对值
  Q_cal <- (1.49/N)*wet.area*hyd.radius^(2/3)*slope^0.5
  return(abs(Q_cal - Q))
}

# 调用optim求解
result <- optim(par = mean(z), 
      fn = Q.calc, 
      lower = min(z), 
      upper = max(z),
      method = 'Brent',
      control = list(maxit = 100000))

运行后得到的result$par就是最优WSE，代入原流量公式可以得到和目标值400接近的计算结果：

# 验证流量准确性
WSE_opt <- result$par
# 代入计算流量
wet.area <- sum(ifelse(pmax(z[-length(z)], z[-1]) > WSE_opt, 0, (x[-1]-x[-length(x)])*(WSE_opt - (z[-1]+z[-length(z)])/2)))
wet.perim <- sum(ifelse(pmax(z[-length(z)], z[-1]) > WSE_opt, 0, sqrt((x[-1]-x[-length(x)])^2 + (z[-1]-z[-length(z)])^2)))
hyd.radius <- wet.area/wet.perim
Q_cal <- (1.49/N)*wet.area*hyd.radius^(2/3)*slope^0.5
print(Q_cal) # 输出值接近400

内容的提问来源于stack exchange，提问作者DFrench208