Why is the objective function in a nonlinear programming (NLP) solver Rsolnp not honored?

Question

The case:

I have a universe of 3 regions (Brazil, New Zealand, USA) (my actual problem is much larger – 31 regions). These three regions are connected through migration. For example, if 10 people move from Brazil to the USA (BRZ-USA) we have immigration to the USA (inflow of people) and emigration from Brazil (outflow of people). I have a data set of migration rates for all possible migration flows in the given universe (3*2 = 6). In addition I have a data set of the population in each region. When I multiply the migration rates by the population, I obtain migrant counts. I can then compute for each region the number of immigrants and the number of emigrants. Subtracting the emigrants from the immigrants results in net migration counts (can be positive or negative). However since we have a balanced system (inflows equal outflows for each region) the sum of the net migrants across all regions should be zero. In addition to the net migrant rates and the population, I also have the net migrant counts from a hypothetical future scenario for each region. But the scenario net migrant counts are different from the counts that I can compute from my data. So I would like to scale the 6 migration rates up and down (by adding or subtracting a fixed number) so that the resulting net migration counts meet the scenario values. I use the nonlinear programming (NLP) solver Rsolnp to accomplish this task (see example script below).

The problem:

I have specified the objective function in the form of a least square equation because it is the goal to force the 6 scalars to be as close to zero as possible. In addition I am using an equality constraint function to meet the scenario values. This all works fine and the solver provides the scalars that I can add to the migration rates leading to the migrant counts to perfectly match the scenario values (see script part “Test if target was reached”). However, I would also like to apply weights (variable: w) to the objective function so that higher values on certain scalars are stronger penalized. However, regardless of how I specify the weights, I always obtain the same solution (see "Example results for different weights"). So it appears that the solver does not honor the objective function. Does anyone have an idea why that is the case and how I could change the objective function so that the use of weights is possible? Thanks a lot for any help!

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

Answer 1

Try different starting values (initScal) than zero; the sum( w * 0^2) = 0 for all w.