带约束的NLS优化：曲线下面积相等的指数拟合问题

我来帮你解决这个带约束的指数拟合问题！首先纠正一个小细节：你之前用log(t)做回归获取初始参数，但指数函数y = a*exp(b*t)的对数形式是log(y) = log(a) + b*t，所以应该用lm(log(y) ~ t)来得到更合理的初始值，这能让优化收敛得更好。

接下来针对你需要的「下行段（t<9）与上行段（t>9）曲线下面积相等」约束，我提供两种可行的实现方案：

我们可以把约束条件和残差平方和目标函数结合成拉格朗日函数，同时优化参数a、b和拉格朗日乘数λ：

y <- c(170, 160, 145, 127, 117, 74, 76, 78, 101, 115, 120, 70, 64, 65)
t <- seq(1,14,1)

# 获取合理初始参数：基于指数函数的对数线性回归
lm_init <- lm(log(y) ~ t)
a_init <- exp(lm_init$coefficients[[1]])
b_init <- lm_init$coefficients[[2]]
lambda_init <- 0  # 拉格朗日乘数初始值设为0

# 定义拉格朗日目标函数
lagrangian_func <- function(pars) {
  a <- pars[1]
  b <- pars[2]
  lambda <- pars[3]
  
  # 计算残差平方和
  fitted <- a * exp(b * t)
  sse <- sum((y - fitted)^2)
  
  # 约束条件：t=1-8与t=10-14的积分面积相等（化简后形式）
  constraint <- exp(14*b) - exp(10*b) - exp(8*b) + exp(b)
  
  # 拉格朗日函数：残差平方和 + 乘数*约束项
  sse + lambda * constraint
}

# 运行优化
result_lagrange <- optim(c(a_init, b_init, lambda_init), lagrangian_func)

# 提取最优参数
a_opt <- result_lagrange$par[1]
b_opt <- result_lagrange$par[2]

# 生成预测值并可视化
pred_opt <- a_opt * exp(b_opt * t)
dat <- data.frame(y=y, t=t, pred=pred_opt)

library(ggplot2)
ggplot(dat, aes(x=t, y=y)) + 
  geom_line(color="black", linewidth=1) + 
  geom_line(aes(y=pred), color="blue", linewidth=1, linetype="dashed") +
  geom_vline(xintercept=9, color="red", linetype="dotted", linewidth=1) +
  labs(title="指数函数拟合（带面积相等约束）", x="t", y="y")

nloptr是专门处理非线性优化的工具包，支持直接设置非线性等式约束，比拉格朗日法更直观：

# 先安装并加载包
install.packages("nloptr")
library(nloptr)

# 目标函数：最小化残差平方和
obj_func <- function(pars) {
  a <- pars[1]
  b <- pars[2]
  fitted <- a * exp(b * t)
  sum((y - fitted)^2)
}

# 定义非线性等式约束：返回值需等于0
eq_constraint <- function(pars) {
  b <- pars[2]
  # 约束的化简形式：两段积分面积相等
  return(exp(14*b) - exp(10*b) - exp(8*b) + exp(b))
}

# 设置优化选项
opts <- list(
  "algorithm" = "NLOPT_LN_COBYLA",  # 适合非线性约束的算法
  "xtol_rel" = 1e-8,
  "maxeval" = 10000
)

# 运行优化
result_nloptr <- nloptr(
  x0 = c(a_init, b_init),
  eval_f = obj_func,
  eval_g_eq = eq_constraint,
  opts = opts
)

# 提取最优参数并生成预测值
a_opt_nl <- result_nloptr$solution[1]
b_opt_nl <- result_nloptr$solution[2]
pred_opt_nl <- a_opt_nl * exp(b_opt_nl * t)

# 可视化
dat_nl <- data.frame(y=y, t=t, pred=pred_opt_nl)
ggplot(dat_nl, aes(x=t, y=y)) + 
  geom_line(color="black", linewidth=1) + 
  geom_line(aes(y=pred), color="green", linewidth=1, linetype="dashed") +
  geom_vline(xintercept=9, color="red", linetype="dotted", linewidth=1) +
  labs(title="指数函数拟合（nloptr带非线性约束）", x="t", y="y")

你可以计算两段的面积来确认约束生效：

# 计算t=1到8的下行段面积
area_down <- (a_opt / b_opt) * (exp(b_opt*8) - exp(b_opt*1))
# 计算t=10到14的上行段面积
area_up <- (a_opt / b_opt) * (exp(b_opt*14) - exp(b_opt*10))

cat("下行段面积：", round(area_down, 2), "\n")
cat("上行段面积：", round(area_up, 2), "\n")
cat("面积差：", round(abs(area_down - area_up), 2), "\n")

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