Shiny应用中solve.QP随机报“约束不一致无解决方案”问题排查

问题背景 #

开发的投资组合优化Shiny应用，运行时约半数概率随机抛出constraints are inconsistent, no solution错误，另一半概率正常。控制台单独执行solve.QP及quadprog相关命令完全正常，仅在Shiny环境中出现问题，怀疑与输入处理或meq约束定义有关，但无法定位具体原因。

核心原因分析 #

1. 全局环境数据污染：getSymbols默认将下载的资产数据存入全局环境，而Shiny的多个输出组件会重复调用dataPrep()响应式函数，导致全局变量被多次覆盖，可能出现数据读取与赋值不同步的情况，进而生成异常的收益率数据、协方差矩阵，引发约束冲突。
2. 协方差矩阵非严格对称：var()生成的协方差矩阵理论上对称，但浮点运算误差可能导致微小不对称，而solve.QP要求输入的Dmat必须严格对称，这种微小差异会随机触发无解错误。
3. 重复代码冗余：QPoptim和meanvarweights函数代码高度重复，增加了出错概率，且不利于维护。

修复方案 #

1. 隔离数据环境：调用getSymbols时指定env = environment()，让数据存入当前响应式的局部环境，避免全局污染和数据冲突。
2. 强制矩阵对称：将协方差矩阵处理为严格对称形式，Dmat = (covmat + t(covmat)) / 2，消除浮点误差导致的不对称问题。
3. 添加错误捕获：在solve.QP调用处添加tryCatch，避免单个优化失败导致整个应用崩溃，同时可输出调试信息。
4. 合并重复函数：将QPoptim和meanvarweights合并为一个函数，同时返回权重和波动率，减少冗余代码。

修改后的完整代码 #

library(quantmod)
library(lubridate)
library(dplyr)
library(data.table)
library(quadprog)
library(shiny)

# Define UI for application that draws a histogram
ui <- fluidPage(
  # Application title
  titlePanel("Robo-Advisor Shiny App"),
  
  # Sidebar with a slider inputs
  fluidRow(
    column(3,
           numericInput(
             inputId = "start",
             label = "Beginning Date (yyymmdd):",
             value = 20160101
           ),
           numericInput(
             inputId = "end",
             label = "Ending Date (yyymmdd):",
             value = 20201231
           ),
           selectInput(
             inputId = "parameter",
             label = "Return optimal portfolio for given:",
             choices = c("mu", "vol"),
             selected = "mu"
           ),
           numericInput(
             inputId = "desired_annual_expected_return",
             label = "Desired annual expected return (in decimal format):",
             value = 0.2
           ), 
           numericInput(
             inputId = "desired_annual_vol",
             label = "Desired annual vol (in decimal format):",
             value = 0.15
           )
    ),
    column(9,align="center",
           fluidRow(
             div(plotOutput("prcPlot")),
             div(tableOutput("titleTable")),
             div(tableOutput("prcTable"))
           )
    )
  ) #close fluidRow
)#close fluidPage

# Define server logic required to draw a histogram
server <- function(input, output) {
  
  #####PREPARING DATA FOR PLOT & TABLE FUNCTIONS
  dataPrep = reactive ({
    #define variables from input boxes--------------------
    startdt = ymd(input$start)
    enddt = ymd(input$end)
    parameter = input$parameter
    d_mean = input$desired_annual_expected_return
    d_sd = input$desired_annual_vol
    
    #download and organize data---------------------------
    symbolList = c("MSFT", "WMT", "AAPL", "IBM", "KO") 
    # 将数据存入当前响应式环境，避免全局污染
    getSymbols(symbolList, from = startdt, to = enddt, src="yahoo", env = environment()) 
    
    #Convert to dataframe
    MSFT = as.data.frame(MSFT) 
    WMT = as.data.frame(WMT)
    AAPL = as.data.frame(AAPL) 
    IBM = as.data.frame(IBM) 
    KO = as.data.frame(KO)
    
    MSFT = to.monthly(MSFT) #converts to monthly frequency
    WMT = to.monthly(WMT) 
    AAPL = to.monthly(AAPL) 
    IBM = to.monthly(IBM)
    KO = to.monthly(KO) 
    
    prices = cbind(MSFT$MSFT.Adjusted, WMT$WMT.Adjusted, AAPL$AAPL.Adjusted, 
                   IBM$IBM.Adjusted, KO$KO.Adjusted)
    
    len = dim(prices)[1]
    returns = as.data.frame(prices[2:len,] / prices[1:(len-1),]) - 1
    names(returns) = c("msft", "wmt", "aapl", "ibm", "ko")
    
    #合并优化函数，同时返回权重和波动率
    meanvarOpt = function(Eport, noshort, N, muvec, covmat){
      ones = array(1,N)
      # 强制Dmat严格对称
      Dmat = (covmat + t(covmat)) / 2
      dvec = array(0,N)
      
      Amat = cbind(muvec, ones)
      b0vec = c(Eport, 1) 
      if(noshort==1) {
        idmat = diag(N)
        Amat = cbind(Amat,idmat)
        b0vec = c(b0vec, array(0,N)) 
      }
      
      # 添加错误捕获，避免单个优化失败导致应用崩溃
      res = tryCatch({
        solve.QP(Dmat, dvec, Amat, b0vec, meq=2)
      }, error = function(e) {
        message("Optimization error: ", e$message)
        return(NULL)
      })
      
      if(is.null(res)){
        return(list(weights = rep(NA, N), vol = NA))
      } else {
        wvec = res$solution
        sigport = sqrt( t(wvec) %*% covmat %*% wvec )
        return(list(weights = wvec, vol = sigport))
      }
    }
    
    #initialize
    N = 5
    muvec = colMeans(returns[,1:N])
    covmat = var(returns[,1:N])
    
    #scroll through Eport values to derive efficient frontier
    mincut = min(muvec)
    maxcut = max(muvec)
    Eportvec = seq(mincut,maxcut,length=300)
    
    # 用lapply调用合并后的函数，提取波动率
    sigportvec = sapply(Eportvec, function(e){
      opt_res = meanvarOpt(e, noshort=1, N, muvec, covmat)
      opt_res$vol
    })
    
    #过滤NA值（如果有优化失败的情况）
    valid_idx = !is.na(sigportvec)
    Eportvec = Eportvec[valid_idx]
    sigportvec = sigportvec[valid_idx]
    
    #annualize stats
    sigportvec = sigportvec * sqrt(12)
    Eportvec = Eportvec * 12
    
    #defining efficient frontier---------------------------------------------------
    if(length(sigportvec) ==0){
      stop("No valid portfolios found in efficient frontier")
    }
    Emincut = Eportvec[which(sigportvec==min(sigportvec))]
    idx = which(Eportvec>=Emincut)
    
    #selecting the optimal portfolio
    w1 = rep(NA, N)
    sig1 = NA
    muportvec = NA
    
    if (parameter=="mu") {
      muportvec = d_mean
      opt_res = meanvarOpt(Eport=muportvec/12, noshort=1, N, muvec, covmat)
      w1 = opt_res$weights
      sig1 = opt_res$vol * sqrt(12)
    } else {
      if(length(idx)==0){
        stop("No valid efficient frontier segment found")
      }
      sigportvec_1 = sigportvec[idx]
      Eportvec_1 = Eportvec[idx]
      a = which(abs(sigportvec_1-d_sd)==min(abs(sigportvec_1-d_sd)))
      muportvec = Eportvec_1[a]
      opt_res = meanvarOpt(Eport=(muportvec/12), noshort=1, N, muvec, covmat)
      w1 = opt_res$weights
      sig1 = opt_res$vol * sqrt(12)
    }
    
    #return data for output function-----------------------
    temp = list(Eportvec = Eportvec, sigportvec = sigportvec, 
                w1=w1, muportvec=muportvec,sig1=sig1,idx=idx)
    
  })
  
  #####Plotting frontier
  output$prcPlot <- renderPlot({
    #calls above function for prepped data------------------------------
    temp = dataPrep()
    sigportvec = temp$sigportvec
    Eportvec = temp$Eportvec
    muportvec = temp$muportvec
    w1 = temp$w1
    sig1 = temp$sig1
    idx= temp$idx
    
    #define variables from input boxes----------------------------------
    startdt = ymd(input$start)
    enddt = ymd(input$end)
    
    #text heading of plot
    startdt_txt = format(startdt, "%Y-%m-%d")
    enddt_txt = format(enddt, "%Y-%m-%d")
    main_text_string = paste("Efficient Frontier (based on data from",startdt_txt,"to",enddt_txt,")")
    
    #plots all minimum variance portfolios------------------------------------------
    plot(sigportvec, Eportvec, 
         xlim=c(0,max(sigportvec, na.rm=T)), ylim=c(0,max(Eportvec, na.rm=T)),
         type="l", xlab="sigma", ylab="E(r)", 
         main=main_text_string, col = "black", lwd=1, lty="dashed")
    
    #now just plot efficient frontier on top------------------------------------------
    if(length(idx)>0){
      lines(x=sigportvec[idx], y=Eportvec[idx], type="l", col = "blue", lwd=2)
    }
    
    #point label text
    if(!is.na(sig1) && !is.na(muportvec)){
      sig1_d = format(round(sig1,2), nsmall = 2)
      muportvec_d = format(round(muportvec,2), nsmall = 2)
      label1 = paste("sig=",sig1_d,", mu=",muportvec_d)
      
      #pick a portfolio along efficient frontier-----------------------
      points(x=c(sig1), y=muportvec, col="blue", lwd=3, pch=1)
      text(x=c(sig1), y=muportvec, labels = c(label1), pos=4)
    }
  })
  
  #Making table heading
  output$titleTable <- renderTable({
    #calls above function for prepped data------------------------------
    temp = dataPrep()
    muportvec = temp$muportvec
    sig1 = temp$sig1
    parameter = input$parameter
    
    #Prepare table heading  format
    if(!is.na(sig1) && !is.na(muportvec)){
      sig1_d = format(round(sig1,2), nsmall = 2)
      muportvec_d = format(round(muportvec,2), nsmall = 2)
      if (parameter=="mu") {
        result = paste("The following portfolio achieves your desired annual mu =", muportvec_d, "with vol =", sig1_d)
      } else {
        result = paste("The following portfolio achieves your desired annual vol =", sig1_d, "with mu =", muportvec_d)
      }
    } else {
      result = "No valid portfolio found for the given parameters"
    }
    
    #making table
    FINALOUT = data.frame(` ` = result)
    FINALOUT
  }, align = 'c')
  
  #####Plotting Summary Table
  output$prcTable <- renderTable({
    #calls above function for prepped data------------------------------
    temp = dataPrep()
    w1 = temp$w1
    
    #reducing decimal places in the weight values
    if(all(!is.na(w1))){
      w1_d = format(round(w1,2),nsmall=2)
      FINALOUT = data.frame(
        Asset = c("MSFT", "WMT", "AAPL", "IBM", "KO"),
        Weight = w1_d
      )
    } else {
      FINALOUT = data.frame(
        Asset = c("MSFT", "WMT", "AAPL", "IBM", "KO"),
        Weight = rep("N/A", 5)
      )
    }
    FINALOUT
  }, align = 'cc')
  
}

# Run the application
shinyApp(ui = ui, server = server)

额外说明 #

• 新增的错误捕获逻辑可以在控制台输出具体的优化错误信息，方便调试。
• 对无效的优化结果（NA值）进行了过滤和处理，避免绘图和表格输出出错。
• 简化了表格生成的代码，去掉了冗余的空行处理，让代码更简洁易读。

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