Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Why does my NLOPT optimization error/fail to solve?

I'm stumped. I have a problem formulated for NLOPT in R. The current problem solves for 180 variables with 28 equality constraints

The code is re-used from a simpler version of the problem, earlier in my script, with 36 variables and 20 equality constraints that solves instantly using NLOPT_LD_SLSQP as the algorithm.

The larger version of the problem with 180 variables produces the following, immediately, when using NLOPT_LD_SLSQP:

NLopt solver status: -4 ( NLOPT_ROUNDOFF_LIMITED: Roundoff errors led 
to a breakdown of the optimization algorithm. In this case, the 
returned minimum may 
still be useful. (e.g. this error occurs in NEWUOA if one tries to 
achieve a tolerance too close to machine precision.) )

This puzzled me, given that it worked on the smaller version of the problem. As well, it returns the starting values and doesn't actually complete any iterations. So I implemented, NLOPT_AUGLAG_LD_EQ as the main algorithm, and NLOPT_LD_SLSQP as the local algorithm. Now the problem fails to solve and produces this:

NLopt solver status: -1 ( NLOPT_FAILURE: Generic failure code. )

If I reduce the tolerance it just fails faster... I took the toy problem into Excel to see if I had failed to formulate properly or if it was somehow infeasible, but it solved immediately. I can give you this file if you want. I took the Excel solution values and populated the functions in R and sure enough my constraints, and objective function appear to be just fine.

I wonder if someone can help me with this problem. Here is some code that will produce the issue (tested and confirmed) for anyone in R:

library(pracma)
library(nloptr)

#My constraint function uses the following:
#RHS of the equality constraints
f.rhs <- c(590.0000,4781.0000,4414.0000,120.0000,224.0000,
849.0000,4693.0000,4374.0000,85.0000,697.0000,0.0000,
0.0000,0.0000,1092.0000,1434.0000,2251.0133,3482.9867,
2316.1813,1873.8187,1622.1450,1206.8550,1240.0000,1233.0000,
933.8532,733.1468,486.7907,395.2093,526.0000)

#matrix of constraints
#the first 10 rows are row total constraints
#the remaining 18 constraints are column total constraints
#this is a 28x180 matrix. It's sort of big to have it in this
#code window, but this code should produce the matrix for you
conmat <- c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.101385354477184,0,0,0,0,0,0,0,0,0,0,0,0,0.101385354477184,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.879602416076347,0,0,0,0,0,0,0,0,0,0,0,0,0,0.879602416076347,0,0,0,0,0,0,0,0,0,0,0,0,0,9.72187831641552,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9.72187831641552,0,0,0,0,0,0,0,0,0,0,0,0,50.825662951981,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,50.825662951981,0,0,0,0,0,0,0,0,0,0,0,61.4161196898944,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,61.4161196898944,0,0,0,0,0,0,0,0,0,0,76.5722969856399,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,76.5722969856399,0,0,0,0,0,0,0,0,0,85.0874667314792,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,85.0874667314792,0,0,0,0,0,0,0,0,70.1228611430807,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,70.1228611430807,0,0,0,0,0,0,0,72.2969630445657,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,72.2969630445657,0,0,0,0,0,0,70.9315070452785,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,70.9315070452785,0,0,0,0,0,54.6520210670868,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,54.6520210670868,0,0,0,0,44.0086494626126,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,44.0086494626126,0,0,0,20.019587567467,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,20.019587567467,0,0,14.1724093345295,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,14.1724093345295,0,10.4922206705268,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,10.4922206705268,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3.56379085643383,0,0,0,0,0,0,0,0,0,0,0,3.56379085643383,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,36.3843642437252,0,0,0,0,0,0,0,0,0,0,0,0,36.3843642437252,0,0,0,0,0,0,0,0,0,0,0,0,0,0,208.581934690648,0,0,0,0,0,0,0,0,0,0,0,0,0,208.581934690648,0,0,0,0,0,0,0,0,0,0,0,0,0,649.993541449925,0,0,0,0,0,0,0,0,0,0,0,0,0,0,649.993541449925,0,0,0,0,0,0,0,0,0,0,0,0,620.425879840303,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,620.425879840303,0,0,0,0,0,0,0,0,0,0,0,532.113517313307,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,532.113517313307,0,0,0,0,0,0,0,0,0,0,487.289086271457,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,487.289086271457,0,0,0,0,0,0,0,0,0,355.571461649492,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,355.571461649492,0,0,0,0,0,0,0,0,370.187737611463,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,370.187737611463,0,0,0,0,0,0,0,377.604110342457,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,377.604110342457,0,0,0,0,0,0,286.391885309974,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,286.391885309974,0,0,0,0,0,230.617639447808,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,230.617639447808,0,0,0,0,150.901768695456,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,150.901768695456,0,0,0,106.827457261503,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,106.827457261503,0,0,102.516716725934,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,102.516716725934,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,478.206023100627,0,0,0,0,0,0,0,0,0,0,478.206023100627,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,656.900445375174,0,0,0,0,0,0,0,0,0,0,0,656.900445375174,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,835.409246572897,0,0,0,0,0,0,0,0,0,0,0,0,835.409246572897,0,0,0,0,0,0,0,0,0,0,0,0,0,0,913.757767618374,0,0,0,0,0,0,0,0,0,0,0,0,0,913.757767618374,0,0,0,0,0,0,0,0,0,0,0,0,0,483.249080637494,0,0,0,0,0,0,0,0,0,0,0,0,0,0,483.249080637494,0,0,0,0,0,0,0,0,0,0,0,0,317.775206732195,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,317.775206732195,0,0,0,0,0,0,0,0,0,0,0,237.647569920133,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,237.647569920133,0,0,0,0,0,0,0,0,0,0,144.341667617821,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,144.341667617821,0,0,0,0,0,0,0,0,0,124.724396017039,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,124.724396017039,0,0,0,0,0,0,0,0,98.7489023828154,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,98.7489023828154,0,0,0,0,0,0,0,55.1727766234815,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,55.1727766234815,0,0,0,0,0,0,44.4279889177619,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,44.4279889177619,0,0,0,0,0,25.9420568421214,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,25.9420568421214,0,0,0,0,18.3650860591977,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,18.3650860591977,0,0,0,18.8065580202466,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,18.8065580202466,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.625528507812058,0,0,0,0,0,0,0,0,0,0.625528507812058,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1.14682340323878,0,0,0,0,0,0,0,0,0,0,1.14682340323878,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,8.19022941749738,0,0,0,0,0,0,0,0,0,0,0,8.19022941749738,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,25.5271393599384,0,0,0,0,0,0,0,0,0,0,0,0,25.5271393599384,0,0,0,0,0,0,0,0,0,0,0,0,0,0,30.1338658434231,0,0,0,0,0,0,0,0,0,0,0,0,0,30.1338658434231,0,0,0,0,0,0,0,0,0,0,0,0,0,33.325484734614,0,0,0,0,0,0,0,0,0,0,0,0,0,0,33.325484734614,0,0,0,0,0,0,0,0,0,0,0,0,34.6159582253122,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,34.6159582253122,0,0,0,0,0,0,0,0,0,0,0,25.9425263509155,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,25.9425263509155,0,0,0,0,0,0,0,0,0,0,28.6721543150432,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,28.6721543150432,0,0,0,0,0,0,0,0,0,24.6001930717219,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,24.6001930717219,0,0,0,0,0,0,0,0,19.4063381820094,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,19.4063381820094,0,0,0,0,0,0,0,15.6269927027328,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,15.6269927027328,0,0,0,0,0,0,7.86354283325046,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,7.86354283325046,0,0,0,0,0,5.56681537403581,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5.56681537403581,0,0,0,0,4.09245697782775,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,4.09245697782775,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.0478232804137659,0,0,0,0,0,0,0,0,0.0478232804137659,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.4097388469824,0,0,0,0,0,0,0,0,0,0.4097388469824,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1.18088960908956,0,0,0,0,0,0,0,0,0,0,1.18088960908956,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2.13050439017778,0,0,0,0,0,0,0,0,0,0,0,2.13050439017778,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2.88537630873373,0,0,0,0,0,0,0,0,0,0,0,0,2.88537630873373,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4.47104432861129,0,0,0,0,0,0,0,0,0,0,0,0,0,4.47104432861129,0,0,0,0,0,0,0,0,0,0,0,0,0,6.7066648766379,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6.7066648766379,0,0,0,0,0,0,0,0,0,0,0,0,7.21292740479352,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,7.21292740479352,0,0,0,0,0,0,0,0,0,0,0,13.2288618280058,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,13.2288618280058,0,0,0,0,0,0,0,0,0,0,17.7096180695658,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,17.7096180695658,0,0,0,0,0,0,0,0,0,32.1794159026695,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,32.1794159026695,0,0,0,0,0,0,0,0,25.9125391288607,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,25.9125391288607,0,0,0,0,0,0,0,39.9666473126856,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,39.9666473126856,0,0,0,0,0,0,28.293474255415,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,28.293474255415,0,0,0,0,0,71.410473393917,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,71.410473393917,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.182608374202251,0,0,0,0,0,0,0,0.182608374202251,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2.83538237187642,0,0,0,0,0,0,0,0,2.83538237187642,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,19.2609656330223,0,0,0,0,0,0,0,0,0,19.2609656330223,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,80.7233037904926,0,0,0,0,0,0,0,0,0,0,80.7233037904926,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,85.6477154102395,0,0,0,0,0,0,0,0,0,0,0,85.6477154102395,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,99.9589781355406,0,0,0,0,0,0,0,0,0,0,0,0,99.9589781355406,0,0,0,0,0,0,0,0,0,0,0,0,0,0,111.246861310206,0,0,0,0,0,0,0,0,0,0,0,0,0,111.246861310206,0,0,0,0,0,0,0,0,0,0,0,0,0,95.4680594823716,0,0,0,0,0,0,0,0,0,0,0,0,0,0,95.4680594823716,0,0,0,0,0,0,0,0,0,0,0,0,110.647972395614,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,110.647972395614,0,0,0,0,0,0,0,0,0,0,0,121.020177240488,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,121.020177240488,0,0,0,0,0,0,0,0,0,0,82.2865998682753,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,82.2865998682753,0,0,0,0,0,0,0,0,0,63.071831445587,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,63.071831445587,0,0,0,0,0,0,0,0,28.5766379506572,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,28.5766379506572,0,0,0,0,0,0,0,26.2025287609293,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,26.2025287609293,0,0,0,0,0,0,24.3533348857096,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,24.3533348857096,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6.77378361061052,0,0,0,0,0,0,6.77378361061052,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,66.8658893692034,0,0,0,0,0,0,0,66.8658893692034,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,301.030835688858,0,0,0,0,0,0,0,0,301.030835688858,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,790.527770127371,0,0,0,0,0,0,0,0,0,790.527770127371,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,571.353574962507,0,0,0,0,0,0,0,0,0,0,571.353574962507,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,465.451417348464,0,0,0,0,0,0,0,0,0,0,0,465.451417348464,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,385.48085698754,0,0,0,0,0,0,0,0,0,0,0,0,385.48085698754,0,0,0,0,0,0,0,0,0,0,0,0,0,0,298.799798756715,0,0,0,0,0,0,0,0,0,0,0,0,0,298.799798756715,0,0,0,0,0,0,0,0,0,0,0,0,0,305.715022425598,0,0,0,0,0,0,0,0,0,0,0,0,0,0,305.715022425598,0,0,0,0,0,0,0,0,0,0,0,0,296.891932828816,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,296.891932828816,0,0,0,0,0,0,0,0,0,0,0,189.45276364351,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,189.45276364351,0,0,0,0,0,0,0,0,0,0,145.213592426377,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,145.213592426377,0,0,0,0,0,0,0,0,0,62.942694262135,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,62.942694262135,0,0,0,0,0,0,0,0,57.7134986817454,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,57.7134986817454,0,0,0,0,0,0,0,35.5644312539887,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,35.5644312539887,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,600.919741246423,0,0,0,0,0,600.919741246423,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,664.169162623053,0,0,0,0,0,0,664.169162623053,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,847.398221128337,0,0,0,0,0,0,0,847.398221128337,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,903.249824396504,0,0,0,0,0,0,0,0,903.249824396504,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,402.366866953328,0,0,0,0,0,0,0,0,0,402.366866953328,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,280.141092267195,0,0,0,0,0,0,0,0,0,0,280.141092267195,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,201.892549678482,0,0,0,0,0,0,0,0,0,0,0,201.892549678482,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,139.225117820461,0,0,0,0,0,0,0,0,0,0,0,0,139.225117820461,0,0,0,0,0,0,0,0,0,0,0,0,0,0,123.454735574923,0,0,0,0,0,0,0,0,0,0,0,0,0,123.454735574923,0,0,0,0,0,0,0,0,0,0,0,0,0,113.363999137925,0,0,0,0,0,0,0,0,0,0,0,0,0,0,113.363999137925,0,0,0,0,0,0,0,0,0,0,0,0,72.6983780508768,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,72.6983780508768,0,0,0,0,0,0,0,0,0,0,0,55.7225581581023,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,55.7225581581023,0,0,0,0,0,0,0,0,0,0,31.2411980593257,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,31.2411980593257,0,0,0,0,0,0,0,0,0,28.6457207488449,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,28.6457207488449,
  0,0,0,0,0,0,0,0,38.4996441218375,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,38.4996441218375,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1.10058560667843,0,0,0,0,1.10058560667843,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3.80905898924575,0,0,0,0,0,3.80905898924575,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,18.1504505832434,0,0,0,0,0,0,18.1504505832434,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,61.5743643896698,0,0,0,0,0,0,0,61.5743643896698,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,52.0421527000771,0,0,0,0,0,0,0,0,52.0421527000771,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,53.6771028451664,0,0,0,0,0,0,0,0,0,53.6771028451664,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,54.7413143789797,0,0,0,0,0,0,0,0,0,0,54.7413143789797,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,46.2299586575942,0,0,0,0,0,0,0,0,0,0,0,46.2299586575942,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,43.5728165091606,0,0,0,0,0,0,0,0,0,0,0,0,43.5728165091606,0,0,0,0,0,0,0,0,0,0,0,0,0,0,36.8659311062071,0,0,0,0,0,0,0,0,0,0,0,0,0,36.8659311062071,0,0,0,0,0,0,0,0,0,0,0,0,0,21.583566050924,0,0,0,0,0,0,0,0,0,0,0,0,0,0,21.583566050924,0,0,0,0,0,0,0,0,0,0,0,0,16.5435811193775,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,16.5435811193775,0,0,0,0,0,0,0,0,0,0,0,6.42202062218569,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6.42202062218569,0,0,0,0,0,0,0,0,0,0,5.88848766417714,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5.88848766417714,0,0,0,0,0,0,0,0,0,3.34484908130525,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3.34484908130525,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.478730062097792,0,0,0,0.478730062097792,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.599532361620929,0,0,0,0,0.599532361620929,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2.08861096083086,0,0,0,0,0,2.08861096083086,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4.67685892604545,0,0,0,0,0,0,4.67685892604545,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6.66062843004317,0,0,0,0,0,0,0,6.66062843004317,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,10.3325985324682,0,0,0,0,0,0,0,0,10.3325985324682,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,17.4366302473174,0,0,0,0,0,0,0,0,0,17.4366302473174,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,23.9406624891688,0,0,0,0,0,0,0,0,0,0,23.9406624891688,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,47.4993402782371,0,0,0,0,0,0,0,0,0,0,0,47.4993402782371,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,75.2636287745429,0,0,0,0,0,0,0,0,0,0,0,0,75.2636287745429,0,0,0,0,0,0,0,0,0,0,0,0,0,0,120.029503824241,0,0,0,0,0,0,0,0,0,0,0,0,0,120.029503824241,0,0,0,0,0,0,0,0,0,0,0,0,0,92.0013786669861,0,0,0,0,0,0,0,0,0,0,0,0,0,0,92.0013786669861,0,0,0,0,0,0,0,0,0,0,0,0,112.914580678886,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,112.914580678886,0,0,0,0,0,0,0,0,0,0,0,103.53378703525,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,103.53378703525,0,0,0,0,0,0,0,0,0,0,216.919314868552,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,216.919314868552)

#Create the matrix from my list of values
conmat <- matrix(conmat,nrow=28,ncol=180)

#Create the constraint function so that it produces 0s
#for the equality constraints

eqn <- function(x){
z=c()
for (i in 1:length(f.rhs)){
  z[i]=x%*%conmat[i,]-f.rhs[i]
}
return(z)
}

#Function for the Jacobian of the constraint function
eqn_grad <- function(x){
jacobian(eqn,x)
}

#Create my objective function: sum of squared error
#Data for the function is in f.obj
f.obj<- c(0,0,0,0.101385354477184,0.879602416076347,9.72187831641552,50.825662951981,61.4161196898944,76.5722969856399,85.0874667314792,70.1228611430807,72.2969630445657,70.9315070452785,54.6520210670868,44.0086494626126,20.019587567467,14.1724093345295,10.4922206705268,0,0,0,3.56379085643383,36.3843642437252,208.581934690648,649.993541449925,620.425879840303,532.113517313307,487.289086271457,355.571461649492,370.187737611463,377.604110342457,286.391885309974,230.617639447808,150.901768695456,106.827457261503,102.516716725934,0,0,0,478.206023100627,656.900445375174,835.409246572897,913.757767618374,483.249080637494,317.775206732195,237.647569920133,144.341667617821,124.724396017039,98.7489023828154,55.1727766234815,44.4279889177619,25.9420568421214,18.3650860591977,18.8065580202466,0,0,0,0.625528507812058,1.14682340323878,8.19022941749738,25.5271393599384,30.1338658434231,33.325484734614,34.6159582253122,25.9425263509155,28.6721543150432,24.6001930717219,19.4063381820094,15.6269927027328,7.86354283325046,5.56681537403581,4.09245697782775,0,0,0,0.0478232804137659,0.4097388469824,1.18088960908956,2.13050439017778,2.88537630873373,4.47104432861129,6.7066648766379,7.21292740479352,13.2288618280058,17.7096180695658,32.1794159026695,25.9125391288607,39.9666473126856,28.293474255415,71.410473393917,0,0,0,0.182608374202251,2.83538237187642,19.2609656330223,80.7233037904926,85.6477154102395,99.9589781355406,111.246861310206,95.4680594823716,110.647972395614,121.020177240488,82.2865998682753,63.071831445587,28.5766379506572,26.2025287609293,24.3533348857096,0,0,0,6.77378361061052,66.8658893692034,301.030835688858,790.527770127371,571.353574962507,465.451417348464,385.48085698754,298.799798756715,305.715022425598,296.891932828816,189.45276364351,145.213592426377,62.942694262135,57.7134986817454,35.5644312539887,0,0,0,600.919741246423,664.169162623053,847.398221128337,903.249824396504,402.366866953328,280.141092267195,201.892549678482,139.225117820461,123.454735574923,113.363999137925,72.6983780508768,55.7225581581023,31.2411980593257,28.6457207488449,38.4996441218375,0,0,0,1.10058560667843,3.80905898924575,18.1504505832434,61.5743643896698,52.0421527000771,53.6771028451664,54.7413143789797,46.2299586575942,43.5728165091606,36.8659311062071,21.583566050924,16.5435811193775,6.42202062218569,5.88848766417714,3.34484908130525,0,0,0,0.478730062097792,0.599532361620929,2.08861096083086,4.67685892604545,6.66062843004317,10.3325985324682,17.4366302473174,23.9406624891688,47.4993402782371,75.2636287745429,120.029503824241,92.0013786669861,112.914580678886,103.53378703525,216.919314868552)

#Objective Function is SSE
fn <- function(x){
sum((f.obj*x-f.obj)^2)
}

#Objective Gradient
fn_grad <- function(x){
grad(fn,x)
}

#Optimization options:
#Starting Values
x0 <- c(matrix(1,1,ncol=length(f.obj)))
#Lower Bound
lb_x <- c(matrix(0,1,ncol=length(f.obj)))
#Upper Bound
ub_x <- c(matrix(2,1,ncol=length(f.obj)))   

Now here is the code to solve with SLSQP:

opts_list <- list("algorithm"="NLOPT_LD_SLSQP", "xtol_rel"=1.0e-8,"maxeval"=1000)

SOLUTION <- nloptr(x0,
                eval_f=fn,
                eval_grad_f=fn_grad,
                lb=lb_x,
                ub=ub_x,
                eval_g_eq = eqn,
                eval_jac_g_eq = eqn_grad,
                opts=opts_list)

SOLUTION

And here is the code to solve with AUGLAG:

local_opts_list <- list("algorithm" = "NLOPT_LD_SLSQP","xtol_rel"=1.0e-8)
opts_list <- list("algorithm"="NLOPT_LD_AUGLAG_EQ", "xtol_rel"=1.0e-8,"maxeval"=1000, "local_opts"=local_opts_list)

SOLUTION <- nloptr(x0,
                eval_f=fn,
                eval_grad_f=fn_grad,
                lb=lb_x,
                ub=ub_x,
                eval_g_eq = eqn,
                eval_jac_g_eq = eqn_grad,
                opts=opts_list)

SOLUTION

Last, but not least, here is a feasible solution produced by Excel:

XL_Sol <- c(1,1,1,0.999840758,0.998619343,0.984736297,0.917305327,0.895766204,0.902054225,0.930449654,0.908095562,0.90555472,0.907109309,0.92775884,0.941019433,0.973119106,0.980882453,1.00818501,1,1,1,0.999277369,0.99265595,0.957137798,1.068940265,0.971466259,0.916751383,1.031659407,1.088943875,1.106986398,1.09298138,1.169120549,1.136193658,1.242584994,1.260102262,1.218743178,1,1,1,1.036812057,1.051986714,1.066969218,0.764543983,0.983061755,1.117743946,1.086522598,1.048114102,1.042021839,1.032920906,1.017781177,1.013537363,1.007869529,1.005459785,1.04562595,1,1,1,0.986367866,0.975008481,0.821522583,0.442429864,0.339768188,0.28380776,0.279905133,0.44151851,0.382930646,0.470415572,0.581960382,0.663070622,0.830426413,0.879918273,0.920434632,1,1,1,0.99976331,0.997972512,0.994156123,0.989331379,0.985345491,0.979153535,0.96871128,0.966148601,0.937972897,0.916916787,0.848764035,0.877734945,0.811485291,0.866305968,0.895221739,1,1,1,0.999728838,0.995792558,0.971414471,0.875797499,0.862210736,0.88101002,0.867466136,0.883358551,0.865529128,0.894845599,0.898520124,0.920974341,0.964087977,0.966910775,1.02107896,1,1,1,1.000481883,1.004812176,0.752446473,1.516261232,1.234618416,1.169698414,1.138623218,1.098571557,1.102118823,1.098211297,1.060360506,1.043658294,1.01881129,1.016899233,1.086075848,1,1,1,0.970757235,0.948975599,1.042336122,0.836515328,0.953263454,1.082324289,1.058181611,1.036019128,1.032395622,1.02937951,1.018031769,1.012837354,1.007153339,1.006385936,1.090482373,1,1,1,0.979339786,0.92850077,0.659360317,0,0.01800471,0.038870008,0,0.145186371,0.194368327,0.318022546,0.600269991,0.693283542,0.880910409,0.890768254,0.94507392,1,1,1,0.998868123,0.998583119,0.995062994,0.988669793,0.983388904,0.978533683,0.963729915,0.94954262,0.90018972,0.84196737,0.764585537,0.80459333,0.760693136,0.77966356,0.896573068)
#eqn(XL_Sol) should give you a vector of numbers that are pretty close to zero.

My question is: why does this model produce the errors I've reported here for each of the SLSQP, and AUGLAG algorithms implemented in the sample code?

Would love some input here. Please let me know if you need any additional information!

like image 526
paperwings Avatar asked Nov 23 '16 18:11

paperwings


1 Answers

An interesting thing about your system is that conmat is rank deficient:

qr(conmat)$rank # 24 < 28 (number of rows in the system)

But the system conmat%*%x = f.rhs is actually consistent and you can find one of its solutions using the answers from this post.

That means that your constraint system has 4 redundant relations. Sometimes issues like this can lead to numerical instability, so let's see what happens if we remove the redundant equations. To do this, we can adapt the answer in this post:

# let's relabel
conmat_raw = conmat
f.rhs_raw = f.rhs

# append rhs as a column to conmat
augmented_matrix <- cbind(conmat, f.rhs)

# qr transform of the transpose
q <- qr(t(augmented_matrix))

# extract a set of linearly dependent rows for the original augmented matrix
reduced <- t(t(augmented_matrix)[,q$pivot[seq(q$rank)]])

# decompose augmented matrix
conmat = reduced[,1:180]
f.rhs = reduced[,181]  

Now, we can run the NLOPT_LD_SLSQP solver. After 1350 iterations, I get a solution. Let's see how it compares to the XL_Sol.

> sol = SOLUTION$solution
> fn(sol)
[1] 58404.02
> fn(XL_Sol)
[1] 317837

So we've improved against the objective. Let's check the solutions using the original constraints:

> sum((conmat_raw%*%sol - f.rhs_raw)^2)
[1] 6.071796e-08
> sum((conmat_raw%*%XL_Sol - f.rhs_raw)^2)
[1] 2.771957e-08
like image 155
Ryan Walker Avatar answered Oct 16 '22 01:10

Ryan Walker