R语言nleqslv函数求解方程组异常问题求助

数据老陈

2026-6-22

R语言nleqslv求解方程组停滞及除数失效问题解决

我在用R语言的nleqslv函数求解10个方程组时，部分输入场景下程序会陷入停滞。通过for循环批量执行求解，引入diviser（除数）优化性能并保证结果合规，但当前示例中除数失效；移除除数后方程组可求解，但结果不正确。完整代码如下：

library(nleqslv)
library(dplyr)

# df_pe1 <-  read.csv("C:\\Users\\IEC\\OneDrive - International Economics Consulting Ltd\\2023\\Partial Equilibrium\\PE test results to WITS SMART\\JG_inputs.csv")

df_pe1 <- data.frame(A= c("A"),
                     B = "B",
                     C = "888",
                     D = "XXX",
                     E = 73792,
                     F = 10.1,
                     G= 1,
                     H = 0.998,
                     I = -7.6004338)

names(df_pe1) <- c("Reporter","Partner","HS6","Commodity","imp_world","imports","initial_tariff","initial_ROW_tariff","imp_elasticity")

df1 <- df_pe1 #%>% filter(HS6 %in% c(80231))

df1 <- df1 %>%
  mutate(diviser = case_when( imports >= imp_world ~ 1*10^(nchar(round(imp_world,0) )+1),
                              imports < imp_world ~ 1*10^(nchar(round(imports,0) )+1)
                              ))

df1 <- df1 %>% 
  mutate(imp_elasticity = case_when(imp_elasticity <= -3 ~ -1.5,
                                    .default = imp_elasticity))

df1[,5] <-df1[,5]/df1$diviser
df1[,6] <-df1[,6]/df1$diviser

df1$ROW_PMxM=df1$imp_world
df1$ROW_TWA=df1$initial_ROW_tariff *100

df1$ROW_Tnew <-df1$initial_ROW_tariff *100

df1$PMxM=df1$imports
df1$TWA<-df1$initial_tariff *100
df1$Tnew <- 0

df1$MELAST =  df1$imp_elasticity #demand elas

df1$XELAST <- 20 #supply elas
df1$PM0<-1

df1$MSUB<- 3#Substitution  els 
#df1$substitution_elasticity

df1$HS<-df1$HS6

df1$M0<-df1$PMxM/df1$PM0
df1$ROW_M0<-df1$ROW_PMxM/df1$PM0

df1$PMA=1

df1$PX0<-df1$PM0/(1+df1$TWA/100)
df1$ROW_PX0<-df1$PM0/(1+df1$ROW_TWA/100)

df1$X0<-df1$M0
df1$ROW_X0<-df1$ROW_M0

df1$MA0<-df1$PMxM+df1$ROW_PMxM
df1$MSUBPAR<--1/df1$MSUB+1

df1$MSHIFT<-df1$MA0/df1$PMA^df1$MELAST

df1$XSHIFT<-df1$X0/df1$PX0^df1$XELAST
df1$ROW_XSHIFT<-df1$ROW_X0/df1$ROW_PX0^df1$XELAST

df1$Mshare<-(1/df1$M0^(df1$MSUBPAR-1))/(1/df1$M0^(df1$MSUBPAR-1)+1/df1$ROW_M0^(df1$MSUBPAR-1))
df1$ROW_Mshare<-(1/df1$ROW_M0^(df1$MSUBPAR-1))/(1/df1$M0^(df1$MSUBPAR-1)+1/df1$ROW_M0^(df1$MSUBPAR-1))

df1$Ashift<-df1$MA0/(df1$Mshare*df1$M0^df1$MSUBPAR+df1$ROW_Mshare*df1$ROW_M0^df1$MSUBPAR)^(1/df1$MSUBPAR)

is.nan.data.frame <- function(x)
  do.call(cbind, lapply(x, is.nan))

df1[is.nan(df1)] <- 0

dftemp<-df1

df<-dftemp

dftemp1<-df

dftemp2<-dftemp1 #1ST BATCH

results<-data.frame(matrix(nrow=nrow(dftemp2), ncol=10))
names(results)<- c("HS","X","PX","ROW_X","ROW_PX","M","ROW_M","PM","ROW_PM","PMA1")
names(dftemp2)
fn3 <- function(x) {
  eq1<-x[1]-df["XSHIFT","Partner"]*x[2]^df["XELAST","Partner"] #X_SUP(P) CIS
  eq2<-x[3]-df["XSHIFT","ROW"]*x[4]^df["XELAST","ROW"] #X_SUP(P) ROW
  eq3<-x[1]-x[5]                                       #EQUIL(P) CIS
  eq4<-x[3]-x[6]                                       #EQUIL(P) ROW
  eq5<-x[7]-x[2]*(1+df["T","Partner"]/100)                 #TARIFF(P) CIS
  eq6<-x[8]-x[4]*(1+df["T","ROW"]/100)                 #TARIFF(P) ROW
  
  eq7<-x[10]-df["ASHIFT",2]*(df["MSHARE","Partner"]*x[5]^df["MSUBPAR",2]+
                               df["MSHARE","ROW"]*x[6]^df["MSUBPAR",2])^(1/df["MSUBPAR",2]) #ARM (Sum ROW and CIS)
  
  eq8<-x[7]-x[9]*df["ASHIFT",2]*(df["MSHARE","Partner"]*(x[5]^df["MSUBPAR",2])+
                                   df["MSHARE","ROW"]*(x[6]^df["MSUBPAR",2]))^((1/df["MSUBPAR",2])-1)*
    df["MSHARE","Partner"]*(x[5]^(df["MSUBPAR",2]-1))                                                      #M_DEM(P) (Sum ROW and CIS) import demand by CIS
  
  eq9<-x[8]-x[9]*df["ASHIFT",2]*(df["MSHARE","Partner"]*(x[5]^df["MSUBPAR",2])+
                                   df["MSHARE","ROW"]*(x[6]^df["MSUBPAR",2]))^((1/df["MSUBPAR",2])-1)*
    df["MSHARE","ROW"]*(x[6]^(df["MSUBPAR",2]-1))                                                      #M_DEM(P) (Sum ROW and CIS) import demand by ROW
  
  eq10<-x[10]- df["MSHIFT",2]*x[9]^df["MELAST",2]    #M_DEMA
  
  return(c(eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10))
}

for (i in 1:nrow(dftemp2)){
  
  print(i)
  # for (i in 4:4){
  df<-data.frame(matrix(nrow=8, ncol=2))
  df[1,1]<-dftemp2[i,"XSHIFT"]
  df[1,2]<-dftemp2[i,"ROW_XSHIFT"]
  df[2,1]<-dftemp2[i,"XELAST"]
  df[2,2]<-dftemp2[i,"XELAST"]
  #df[3,1]<-dftemp2[i,"TWA"]
  df[3,1]<-dftemp2[i,"Tnew"]
  #df[3,2]<-dftemp2[i,"W_TWA"]
  df[3,2]<-dftemp2[i,"ROW_Tnew"]
  df[4,1]<-dftemp2[i,"Ashift"]
  df[4,2]<-dftemp2[i,"Ashift"]
  df[5,1]<-dftemp2[i,"Mshare"]
  df[5,2]<-dftemp2[i,"ROW_Mshare"]
  
  df[6,1]<-dftemp2[i,"MSUBPAR"]
  df[6,2]<-dftemp2[i,"MSUBPAR"]
  
  df[7,1]<-dftemp2[i,"MSHIFT"]
  df[7,2]<-dftemp2[i,"MSHIFT"]
  
  df[8,1]<-dftemp2[i,"MELAST"]
  df[8,2]<-dftemp2[i,"MELAST"]
  
  rownames(df)<-c("XSHIFT","XELAST","T","ASHIFT","MSHARE","MSUBPAR","MSHIFT","MELAST")
  names(df)<-c("Partner","ROW")
  
  xstart<-c(1,1,1,1,1,1,1,1,1,1)
  
  # mxit <- 10^8*dftemp2[1,"diviser"]
  z <- nleqslv(xstart, fn3, control=list(maxit=10^8,allowSingular=TRUE)
  )
  
  results[i,"HS"]<-dftemp2[i,"HS"]
  results[i,"X"]<-z$x[[1]]
  results[i,"PX"]<-z$x[[2]]
  results[i,"ROW_X"]<-z$x[[3]]
  results[i,"ROW_PX"]<-z$x[[4]]
  
  results[i,"M"]<-z$x[[5]]
  results[i,"ROW_M"]<-z$x[[6]]
  results[i,"PM"]<-z$x[[7]]
  results[i,"ROW_PM"]<-z$x[[8]]
  results[i,"PMA1"]<-z$x[[9]]
  
}

df<-full_join(dftemp2,results,by=c("HS"))

df$change_XS<-(df$XSHIFT/(df$XELAST+1)*(df$PX^(df$XELAST+1)-df$PX0^(df$XELAST+1))) * df$diviser #CH_XGFT(P) CIS
df$ROW_change_XS<-(df$ROW_XSHIFT/(df$XELAST+1)*(df$ROW_PX^(df$XELAST+1)-df$ROW_PX0^(df$XELAST+1))) * df$diviser #CH_XGFT(P) ROW

df$MS_change<- df$diviser * (ifelse((df$MELAST+1) != 0, df$MSHIFT/(df$MELAST+1), 0)*(df$PMA^(df$MELAST+1)-df$PMA1^(df$MELAST+1)))   #CH_MGFT #Welfare effect

tnew <- df1$Tnew[[1]]

df$MELAST <- df$MELAST

#Initial Imports
df$imports <- df$imports *  df$diviser 

df$imp_world <- df$imp_world *  df$diviser  #current import ROW

#New imports
df$M <- df$M * df$diviser #new import 

df$ROW_M <- df$ROW_M * df$diviser #new import ROW

df$change_import <- df$M - df$imports #change in import from country 

df$ROW_change_import <- df$ROW_M - df$imp_world # change in import ROW

#tariff revenue
df$TR_change <- ((df$imports * (df$Tnew/100 - df$TWA/100 )) + 
                   (df$ROW_change_import * df$ROW_TWA/100)) 

df$TR0 <- df$TWA/100 * (df$imports + df$imp_world)

df$TR1 <-  df$TR0 + df$TR_change

问题根源分析

除数缩放逻辑不完整：仅对imp_world和imports做了除法缩放，但未同步调整依赖这些值计算的中间变量（如M0、MA0、XSHIFT等），导致方程组变量尺度不匹配，求解器无法收敛。
初始值与变量尺度不匹配：缩放后变量值极小，但初始值xstart仍设为全1，与实际变量尺度差距过大，nleqslv迭代时难以找到收敛方向，陷入停滞。
移除除数后的数值精度问题：原始变量值过大，幂运算（如x[2]^20）产生超出精度范围的极值，导致求解结果失真。

针对性解决方案

1. 统一变量缩放逻辑

所有依赖imp_world和imports的中间变量都按diviser缩放，确保方程组内变量尺度一致：

df1 <- df1 %>%
  mutate(
    imp_world = imp_world / diviser,
    imports = imports / diviser,
    PMxM = imports,
    ROW_PMxM = imp_world,
    M0 = PMxM / PM0,
    ROW_M0 = ROW_PMxM / PM0,
    MA0 = PMxM + ROW_PMxM
  )

2. 调整初始值匹配缩放尺度

基于缩放后的变量值设置合理初始值，缩小与真实解的差距：

xstart <- c(
  dftemp2[i,"XSHIFT"],
  dftemp2[i,"PX0"],
  dftemp2[i,"ROW_XSHIFT"],
  dftemp2[i,"ROW_PX0"],
  dftemp2[i,"M0"],
  dftemp2[i,"ROW_M0"],
  dftemp2[i,"PM0"]*(1+dftemp2[i,"Tnew"]/100),
  dftemp2[i,"PM0"]*(1+dftemp2[i,"ROW_Tnew"]/100),
  dftemp2[i,"PMA"],
  dftemp2[i,"MSHIFT"]
)

3. 优化方程组函数的参数传递

避免循环中df被覆盖的歧义，改用参数列表传递当前行数据：

fn3 <- function(x, params) {
  eq1 <- x[1] - params$XSHIFT * x[2]^params$XELAST
  eq2 <- x[3] - params$ROW_XSHIFT * x[4]^params$XELAST
  eq3 <- x[1] - x[5]
  eq4 <- x[3] - x[6]
  eq5 <- x[7] - x[2]*(1+params$Tnew/100)
  eq6 <- x[8] - x[4]*(1+params$ROW_Tnew/100)
  
  term <- params$Mshare*x[5]^params$MSUBPAR + params$ROW_Mshare*x[6]^params$MSUBPAR
  eq7 <- x[10] - params$Ashift*(term)^(1/params$MSUBPAR)
  
  eq8 <- x[7] - x[9]*params$Ashift*(term)^((1/params$MSUBPAR)-1)*params$Mshare*(x[5]^(params$MSUBPAR-1))
  
  eq9 <- x[8] - x[9]*params$Ashift*(term)^((1/params$MSUBPAR)-1)*params$ROW_Mshare*(x[6]^(params$MSUBPAR-1))
  
  eq10 <- x[10] - params$MSHIFT*x[9]^params$MELAST
  
  return(c(eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10))
}

4. 调整求解器控制参数

降低maxit至合理值，避免无效迭代；尝试更换求解算法提升收敛效率：

z <- nleqslv(xstart, fn3, params=params, control=list(maxit=1000, allowSingular=TRUE, method="Newton"))

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

R语言nleqslv函数求解方程组异常问题求助

问题根源分析 #

针对性解决方案 #

1. 统一变量缩放逻辑 #

2. 调整初始值匹配缩放尺度 #

3. 优化方程组函数的参数传递 #