0

案例:爲什麼非線性規劃(NLP)解算器Rsolnp中的目標函數不符合?

我有3個地區(巴西,新西蘭,美國)的宇宙(我的實際問題是更大 - 31個區)。這三個地區通過遷移相連。例如,如果有10人從巴西轉移到美國(BRZ-USA),我們就有移民到美國(人流入)和從巴西移民(人流出)。我有一個關於給定宇宙中所有可能的遷移流的遷移率數據集(3 * 2 = 6)。另外,我還有每個地區的人口數據集。當我將遷移率乘以人口時,我獲得了移民數量。然後我可以計算每個地區的移民人數和移民人數。從移民中扣除移民導致淨移民數量(可以是正數或負數)。然而,由於我們有一個平衡的制度(每個地區的流入流量相等),所有地區的淨流動人口總和應爲零。除了淨移民率和人口外,我還會從每個地區假設的未來情景中獲得淨移民數量。但情景淨移民數量與我可以從我的數據計算的數量不同。因此,我想通過增加或減少一個固定數字來擴大和縮小6個遷移率,以使得到的淨遷移數量符合方案值。我使用非線性規劃(NLP)解算器Rsolnp來完成此任務(請參閱下面的示例腳本)。

問題:

我已指定在最小二乘方程的形式的目標函數,因爲它是迫使6個定標器以儘可能接近零越好目標。另外,我正在使用等式約束函數來滿足場景值。這一切都正常工作,解決方案提供了可以添加到遷移率的標量,從而導致遷移計數與場景值完美匹配(請參閱腳本部分「測試目標是否已達到」)。但是,我還想將權重(變量:w)應用於目標函數,以便某些標量的較高值受到較大的懲罰。但是,無論我如何指定權重,我總是會獲得相同的解決方案(請參閱「不同權重的示例結果」)。所以看起來求解器不遵守目標函數。有沒有人有一個想法,爲什麼是這種情況,以及如何改變目標函數,以便使用權重是可能的?非常感謝您的幫助!

library(Rsolnp) 

# Regions 
regUAll=c("BRZ","NZL","USA") # "BRZ"=Brazil; "NZL"=New Zealand; "USA"=United States 

#* Generate unique combinations of regions 
uCombi=expand.grid(regUAll,regUAll,stringsAsFactors=F) 
uCombi=uCombi[uCombi$Var1!=uCombi$Var2,] # remove same region combination (e.g., BRZ-BRZ) 
uCombi=paste(uCombi$Var2,uCombi$Var1,sep="-") 

#* Generate data frames 
# Migration rates - rows represent major age groups (row1=0-25 years, row2=26-50 years, row3=51-75 years) 
dfnm=data.frame(matrix(rep(c(0.01,0.04,0.02),length(uCombi)),ncol=length(uCombi),nrow=3),stringsAsFactors=F) # generate empty df 
names(dfnm)=uCombi # assign variable names 

# Population (number of people) in region of origin 
pop=c(rep(c(20,40,10),2),rep(c(4,7,2),2),rep(c(30,70,50),2)) 
dfpop=data.frame(matrix(pop,ncol=length(uCombi),nrow=3),stringsAsFactors=F) # generate empty df 
names(dfpop)=uCombi # assign variable names 

#* Objective function for optimization 
# Note: Least squares method to keep the additive scalers as close to 0 as possible 
#  The sum expression allows for flexible numbers of scalars to be included but is identical to: w[1](scal[1]-0)^2+w[2](scal[2]-0)^2+w[3](scal[3]-0)^2+w[4](scal[4]-0)^2+w[5](scal[5]-0)^2+w[6](scal[6]-0)^2 
f.main=function(scal,nScal,w,dfnm,dfpop,regUAll){ 
    sum(w*(scal[1:nScal]-0)^2) 
} 

#* Equality contraint function 
f.equal=function(scal,nScal,w,dfnm,dfpop,regUAll){ 

    #* Adjust net migration rates by scalar 
    for(s in 1:nScal){ 
    dfnm[,s]=dfnm[,s]+scal[s] 
    } 

    #* Compute migration population from data 
    nmp=sapply(dfpop*dfnm,sum) # sums migration population across age groups 

    nmd=numeric(length(regUAll)); names(nmd)=regUAll      # generate named vector to be filled with values 
    for(i in 1:length(regUAll)){ 
    colnEm=names(nmp)[grep(paste0("^",regUAll[i],"-.*"),names(nmp))]  # emigration columns 
    colnIm=names(nmp)[grep(paste0("^.*","-",regUAll[i],"$"),names(nmp))] # immigration columns 
    nmd[regUAll[i]]=sum(nmp[colnIm])-sum(nmp[colnEm])      # compute net migration population = immigration - emigration 
    } 
    nmd=nmd[1:(length(nmd)-1)] # remove the last equality constraint value - not needed because we have a closed system in which global net migration=0 

    return(nmd) 
} 

#* Set optimization parameters 
cpar2=list(delta=1,tol=1,outer.iter=10,trace=1) # optimizer settings 
nScal=ncol(dfnm)     # number of scalars to be used 
initScal=rep(0,nScal)    # initial values of additive scalars 
lowScal=rep(-1,nScal)    # lower bounds on scalars 
highScal=rep(1,nScal)    # upper bounds on scalars 
nms=c(-50,10)      # target values: BRZ=-50, NZL=10, USA=40; last target value does not need to be included since we deal with a closed system in which global net migration sums to 0 
w=c(1,1,1,1,1,1) # unity weights 
#w=c(1,1,2,2,1,1) # double weight on NZL 
#w=c(5,1,2,7,1,0.5) # mixed weights 


#* Perform optimization using solnp 
solRes=solnp(initScal,fun=f.main,eqfun=f.equal,eqB=nms,LB=lowScal,UB=highScal,control=cpar2, 
      nScal=nScal,w=w,dfnm=dfnm,dfpop=dfpop,regUAll=regUAll) 

scalSol=solRes$pars # return optimized values of scalars 

# Example results for different weights 
#[1] 0.101645349 0.110108019 -0.018876993 0.001571639 -0.235945755 -0.018134294 # w=c(1,1,1,1,1,1) 
#[1] 0.101645349 0.110108019 -0.018876993 0.001571639 -0.235945755 -0.018134294 # w=c(1,1,2,2,1,1) 
#[1] 0.101645349 0.110108019 -0.018876993 0.001571639 -0.235945755 -0.018134294 # w=c(5,1,2,7,1,0.5) 

#*** Test if target was reached 
# Adjust net migration rates using the optimized scalars 
for(s in 1:nScal){ 
    dfnm[,s]=dfnm[,s]+scalSol[s] 
} 

# Compute new migration population 
nmp=sapply(dfpop*dfnm,sum) # sums migration population across age groups 

nmd=numeric(length(regUAll)); names(nmd)=regUAll      # generate named vector to be filled with values 
for(i in 1:length(regUAll)){ 
    colnEm=names(nmp)[grep(paste0("^",regUAll[i],"-.*"),names(nmp))]  # emigration columns 
    colnIm=names(nmp)[grep(paste0("^.*","-",regUAll[i],"$"),names(nmp))] # immigration columns 
    nmd[regUAll[i]]=sum(nmp[colnIm])-sum(nmp[colnEm])      # compute net migration population = immigration - emigration 
} 

nmd # should be -50,10,40 if scalars work correctly 
+0

嘗試不同於零的起始值('initScal');對於所有w,和(w * 0^2)= 0。 – 2014-12-01 19:10:10

+0

偉大的建議@NealFultz!當將初始值設置爲不同於零的值時,權重實際上就起作用。如果您將答案作爲普通帖子發佈,我可以接受。唯一令我擔心的是初始值的值對所得到的標量值有很大影響。 – Raphael 2014-12-02 16:18:16

回答

0

嘗試不同的起始值(initScal)大於零;所有wsum(w * 0^2) = 0

+0

對於初始值的小數值的選擇,重要的是要考慮在選項(cpar)中設置的步長(delta)和容差(tol)。當使用小的起始值和相對粗略的步長和容差時(例如,起始值0.0001;增量= 0.1; tol = 0.1),目標函數和權重不受尊重。但是,當增加步長和容差的靈敏度時,可以使用較小的起始值(例如,起始值0.0001; delta = 0.001; tol = 0.001),並且目標函數和權重仍然有效。 – Raphael 2014-12-02 18:50:06