Rajarshi Bhadra
Reputation: 1944
Using Rsolnp for optimization
I have data frame which looks like this
> Data_bounded
Key channel oas spend beta alpha LB UB Lower Upper
1 DRM_BBY DRM BBY 3.599127e+06 0.4355527 3.790177e+04 30 30 2.519389e+06 4.678865e+06
2 DRM_ESS DRM ESS 1.560708e+07 1.0000000 0.000000e+00 30 30 1.092496e+07 2.028921e+07
3 DRM_STG DRM STG 2.053543e+06 0.4846562 6.602477e+04 30 30 1.437480e+06 2.669605e+06
4 DRM_UNU DRM UNU 2.541332e+06 0.5000000 1.296799e+03 30 30 1.778932e+06 3.303732e+06
5 EML_BBY EML BBY 0.000000e+00 1.0000000 0.000000e+00 0 0 0.000000e+00 0.000000e+00
6 EML_BTY EML BTY 0.000000e+00 1.0000000 0.000000e+00 0 0 0.000000e+00 0.000000e+00
7 EML_SHL EML SHL 0.000000e+00 1.0000000 0.000000e+00 0 0 0.000000e+00 0.000000e+00
8 ETC_BTY ETC BTY 0.000000e+00 1.0000000 0.000000e+00 100 10 0.000000e+00 0.000000e+00
9 ETC_ESS ETC ESS 3.040605e+06 0.5430741 3.661138e+03 100 10 0.000000e+00 3.344666e+06
10 ETC_LAM ETC LAM 2.184844e+07 0.5570489 2.379726e+04 100 10 0.000000e+00 2.403328e+07
11 ETN_ESS ETN ESS 2.571237e+06 0.5122024 1.400936e+03 100 10 0.000000e+00 2.828361e+06
12 ETN_LAM ETN LAM 3.446503e+07 0.6619932 5.246261e+03 100 10 0.000000e+00 3.791154e+07
13 ETS_BTY ETS BTY 0.000000e+00 1.0000000 0.000000e+00 100 10 0.000000e+00 0.000000e+00
14 ETS_ESS ETS ESS 3.179937e+05 0.5747293 4.705728e+03 100 10 0.000000e+00 3.497930e+05
15 ETS_LAM ETS LAM 4.071746e+06 0.7597932 2.232724e+03 100 10 0.000000e+00 4.478921e+06
16 FSI_BTY FSI BTY 7.093370e+04 1.0000000 0.000000e+00 5 60 6.738701e+04 1.134939e+05
17 FSI_STG FSI STG 5.852890e+05 0.3000000 1.216848e+06 5 60 5.560246e+05 9.364624e+05
18 MAG_BBY MAG BBY 1.591777e+06 0.5816724 1.753382e+04 10 5 1.432599e+06 1.671366e+06
19 MAG_BTY MAG BTY 8.526921e+06 0.3000000 1.704103e+05 10 5 7.674229e+06 8.953267e+06
20 MAG_SHL MAG SHL 5.546460e+05 0.6597280 1.604735e+03 10 5 4.991814e+05 5.823783e+05
21 MAG_UNU MAG UNU 6.349920e+05 0.7000000 2.145462e+02 10 5 5.714928e+05 6.667416e+05
22 MSM_BTY MSM BTY 5.953225e+06 1.0000000 0.000000e+00 30 30 4.167258e+06 7.739193e+06
23 MSM_CLI MSM CLI 3.095554e+05 1.0000000 0.000000e+00 30 30 2.166888e+05 4.024220e+05
24 MSM_STG MSM STG 2.051692e+06 0.5000000 1.161424e+03 30 30 1.436184e+06 2.667199e+06
25 OND_BBY OND BBY 1.195467e+06 0.4000000 1.279868e+03 20 20 9.563732e+05 1.434560e+06
26 OND_BTY OND BTY 7.513692e+05 0.3048512 3.194817e+05 20 20 6.010953e+05 9.016430e+05
27 OND_CLI OND CLI 1.901847e+05 0.5009272 1.597470e+03 20 20 1.521478e+05 2.282217e+05
28 OND_ESS OND ESS 8.269617e+05 0.3000000 1.479966e+05 20 20 6.615694e+05 9.923540e+05
29 OND_SHL OND SHL 2.844579e+06 0.5606207 5.739453e+02 20 20 2.275663e+06 3.413495e+06
30 OND_STG OND STG 6.673592e+04 1.0000000 0.000000e+00 20 20 5.338874e+04 8.008310e+04
31 OND_UNU OND UNU 5.940489e+05 0.5097571 1.905513e+03 20 20 4.752391e+05 7.128586e+05
32 OOH_BBY OOH BBY 3.149710e+05 1.0000000 0.000000e+00 10 5 2.834739e+05 3.307195e+05
33 OOH_BTY OOH BTY 5.000000e-07 1.0000000 0.000000e+00 10 5 4.500000e-07 5.250000e-07
34 OOH_CLI OOH CLI 1.299990e+05 0.3000000 1.657615e+04 10 5 1.169991e+05 1.364990e+05
35 OOH_STG OOH STG 7.059643e+02 1.0000000 0.000000e+00 10 5 6.353679e+02 7.412625e+02
36 PSS_BBY PSS BBY 9.281210e+05 0.7092288 2.608441e+02 10 10 8.353089e+05 1.020933e+06
37 PSS_BTY PSS BTY 7.408000e+04 0.4885998 1.246914e+04 10 10 6.667200e+04 8.148800e+04
38 PSS_CLI PSS CLI 7.849900e+04 0.7897747 1.361035e+02 10 10 7.064910e+04 8.634890e+04
39 PSS_SHL PSS SHL 1.003498e+06 0.5183048 8.555700e+03 10 10 9.031482e+05 1.103848e+06
40 PSS_STG PSS STG 3.353000e+05 0.3000000 3.811032e+05 10 10 3.017700e+05 3.688300e+05
41 PSS_UNU PSS UNU 2.394850e+05 0.4749274 1.651913e+03 10 10 2.155365e+05 2.634335e+05
42 RAD_BBY RAD BBY 9.182681e+05 0.5000000 1.979574e+03 10 10 8.264413e+05 1.010095e+06
43 RAD_ESS RAD ESS 0.000000e+00 1.0000000 0.000000e+00 10 10 0.000000e+00 0.000000e+00
44 RAD_STG RAD STG 2.852897e+06 0.4176579 4.354847e+04 10 10 2.567607e+06 3.138186e+06
45 RCM_BTY RCM BTY 1.246777e+05 0.4000000 7.491181e+03 0 15 1.246777e+05 1.433794e+05
46 RCM_ESS RCM ESS 2.242704e+04 1.0000000 0.000000e+00 0 15 2.242704e+04 2.579110e+04
47 RCM_FOB RCM FOB 1.094317e+06 0.3684433 1.062531e+05 0 15 1.094317e+06 1.258465e+06
48 ROP_CLI ROP CLI 1.564082e+05 1.0000000 0.000000e+00 20 5 1.251266e+05 1.642286e+05
49 ROP_SHL ROP SHL 1.512256e+05 0.6284284 2.321274e+01 20 5 1.209804e+05 1.587868e+05
50 ROP_STG ROP STG 6.431417e+06 0.4100540 1.510802e+04 20 5 5.145133e+06 6.752987e+06
51 STC_ESS STC ESS 0.000000e+00 1.0000000 0.000000e+00 100 10 0.000000e+00 0.000000e+00
52 STC_LAM STC LAM 2.822590e+05 0.6630324 6.334003e+03 100 10 0.000000e+00 3.104849e+05
53 STE_LAM STE LAM 1.023410e+06 0.6736357 4.057302e+03 10 5 9.210691e+05 1.074581e+06
54 STN_ESS STN ESS 0.000000e+00 1.0000000 0.000000e+00 100 10 0.000000e+00 0.000000e+00
55 STN_LAM STN LAM 4.273955e+06 0.7610039 3.018853e+03 100 10 0.000000e+00 4.701351e+06
56 STS_ESS STS ESS 0.000000e+00 1.0000000 0.000000e+00 10 5 0.000000e+00 0.000000e+00
57 STS_LAM STS LAM 5.558187e+05 0.4623422 5.090183e+04 10 5 5.002369e+05 5.836097e+05
58 VDO_ESS VDO ESS 2.206729e+06 0.8000000 2.697051e+02 0 15 2.206729e+06 2.537739e+06
59 VDO_LAM VDO LAM 2.759217e+06 0.7601210 8.131849e+01 0 15 2.759217e+06 3.173100e+06
60 VDO_UNU VDO UNU 0.000000e+00 1.0000000 0.000000e+00 0 15 0.000000e+00 0.000000e+00
I am using the optimization package Rsolnp
Now my objective is to minimize the function:
-sum(alpha*Spend^beta)
subject to the constraint spend lies between Lower and Upper and that sum of optimum spend values be equal to sum of column spend. However the solution is not yielding any result.Any help will be greatly appreciated. The code I am using is
Data_bounded is my data.
alpha<-Data_bounded$alpha
beta<-Data_bounded$beta
LB_data<-Data_bounded$Lower
UB_data<-Data_bounded$Upper
total_spend<-sum(Data_bounded$spend)
fn1<-function(x)
{
z<-seq(0,0,length=nrow(Data_bounded))
for(i in (1:nrow(Data_bounded))){
z[i]<- alpha[i]*(x[i])^beta[i]
}
return(-sum(z))
}
eq1<-function(x)
{
sum(x)
}
x0 = Data_bounded$spend
opti<-solnp(x0,fun=fn1,LB=LB_data,UB=UB_data,eqfun=eq1,eqB=total_spend)
optimum_spend<-opti$pars
> SAS_soln
channel oas spend beta alpha key UpperBd LowerBd Optimum_Spend opt_Inc_Sales
1 DRM BBY 3599126.68 0.43555 37901.77 DRM_BBY 4678864.68 2519388.67 4678864.68 30468915.4
2 DRM ESS 15607084.43 1.00000 0.00 DRM_ESS 20289209.76 10924959.10 10924959.10 0
3 DRM STG 2053542.67 0.48466 66024.77 DRM_STG 2669605.47 1437479.87 2669605.47 85965546.97
4 DRM UNU 2541332.00 0.50000 1296.80 DRM_UNU 3303731.59 1778932.40 1778932.40 1729627.27
5 EML BBY 0.00 1.00000 0.00 EML_BBY 0.00 0.00 0.00 .
6 EML BTY 0.00 1.00000 0.00 EML_BTY 0.00 0.00 0.00 .
7 EML SHL 0.00 1.00000 0.00 EML_SHL 0.00 0.00 0.00 .
8 ETC BTY 0.00 1.00000 0.00 ETC_BTY 0.00 0.00 0.00 .
9 ETC ESS 3040605.45 0.54307 3661.14 ETC_ESS 3344666.00 0.00 3344666.00 12788595.53
10 ETC LAM 21848435.35 0.55705 23797.26 ETC_LAM 24033278.89 0.00 24033278.89 307610216.1
11 ETN ESS 2571237.02 0.51220 1400.94 ETN_ESS 2828360.72 0.00 1436429.44 1996155.78
12 ETN LAM 34465032.77 0.66199 5246.26 ETN_LAM 37911536.05 0.00 37911536.05 545688932
13 ETS BTY 0.00 1.00000 0.00 ETS_BTY 0.00 0.00 0.00 .
14 ETS ESS 317993.68 0.57473 4705.73 ETS_ESS 349793.05 0.00 349793.05 7224665.67
15 ETS LAM 4071746.49 0.75979 2232.72 ETS_LAM 4478921.14 0.00 4478921.14 252552248.9
16 FSI BTY 70933.70 1.00000 0.00 FSI_BTY 113493.92 67387.02 67387.02 0
17 FSI STG 585289.02 0.30000 1216848.42 FSI_STG 936462.43 556024.57 936462.43 75280689.89
18 MAG BBY 1591777.00 0.58167 17533.82 MAG_BBY 1671365.85 1432599.30 1671365.85 73057646.57
19 MAG BTY 8526921.00 0.30000 170410.33 MAG_BTY 8953267.05 7674228.90 8660097.95 20547204.92
20 MAG SHL 554646.00 0.65973 1604.74 MAG_SHL 582378.30 499181.40 582378.30 10206360.44
21 MAG UNU 634992.00 0.70000 214.55 MAG_UNU 666741.60 571492.80 666741.60 2560301.93
22 MSM BTY 5953225.18 1.00000 0.00 MSM_BTY 7739192.73 4167257.63 4167257.63 0
23 MSM CLI 309555.42 1.00000 0.00 MSM_CLI 402422.05 216688.80 216688.80 0
24 MSM STG 2051691.88 0.50000 1161.42 MSM_STG 2667199.44 1436184.32 1436184.32 1391860.93
25 OND BBY 1195466.55 0.40000 1279.87 OND_BBY 1434559.86 956373.24 956373.24 315802.98
26 OND BTY 751369.17 0.30485 319481.69 OND_BTY 901643.00 601095.33 901643.00 20885532.25
27 OND CLI 190184.72 0.50093 1597.47 OND_CLI 228221.67 152147.78 228221.67 771932.76
28 OND ESS 826961.70 0.30000 147996.58 OND_ESS 992354.04 661569.36 992354.04 9316476.3
29 OND SHL 2844578.87 0.56062 573.95 OND_SHL 3413494.65 2275663.10 2275663.10 2102793.14
30 OND STG 66735.92 1.00000 0.00 OND_STG 80083.10 53388.74 53388.74 0
31 OND UNU 594048.85 0.50976 1905.51 OND_UNU 712858.62 475239.08 712858.62 1834941.14
32 OOH BBY 314971.00 1.00000 0.00 OOH_BBY 330719.55 283473.90 283473.90 0
33 OOH BTY 0.00 1.00000 0.00 OOH_BTY 0.00 0.00 0.00 .
34 OOH CLI 129999.00 0.30000 16576.15 OOH_CLI 136498.95 116999.10 136498.95 575469.83
35 OOH STG 705.96 1.00000 0.00 OOH_STG 741.26 635.37 635.37 0
36 PSS BBY 928121.00 0.70923 260.84 PSS_BBY 1020933.10 835308.90 1020933.10 4765788.27
37 PSS BTY 74080.00 0.48860 12469.14 PSS_BTY 81488.00 66672.00 81488.00 3128929.49
38 PSS CLI 78499.00 0.78977 136.10 PSS_CLI 86348.90 70649.10 86348.90 1077457.67
39 PSS SHL 1003498.00 0.51830 8555.70 PSS_SHL 1103847.80 903148.20 1103847.80 11596424.97
40 PSS STG 335300.00 0.30000 381103.19 PSS_STG 368830.00 301770.00 368830.00 17827500.41
41 PSS UNU 239485.00 0.47493 1651.91 PSS_UNU 263433.50 215536.50 263433.50 620029.25
42 RAD BBY 918268.08 0.50000 1979.57 RAD_BBY 1010094.89 826441.27 1010094.89 1989541.09
43 RAD ESS 0.00 1.00000 0.00 RAD_ESS 0.00 0.00 0.00 .
44 RAD STG 2852896.79 0.41766 43548.47 RAD_STG 3138186.47 2567607.11 3138186.47 22509301
45 RCM BTY 124677.72 0.40000 7491.18 RCM_BTY 143379.38 124677.72 143379.38 865256.72
46 RCM ESS 22427.04 1.00000 0.00 RCM_ESS 25791.10 22427.04 22427.04 0
47 RCM FOB 1094317.32 0.36844 106253.06 RCM_FOB 1258464.92 1094317.32 1258464.92 18783966.92
48 ROP CLI 156408.20 1.00000 0.00 ROP_CLI 164228.61 125126.56 125126.56 0
49 ROP SHL 151225.56 0.62843 23.21 ROP_SHL 158786.84 120980.45 120980.45 36295.19
50 ROP STG 6431416.66 0.41005 15108.02 ROP_STG 6752987.49 5145133.33 5145133.33 8535764.15
51 STC ESS 0.00 1.00000 0.00 STC_ESS 0.00 0.00 0.00 .
52 STC LAM 282258.99 0.66303 6334.00 STC_LAM 310484.89 0.00 310484.89 27738245.2
53 STE LAM 1023410.08 0.67364 4057.30 STE_LAM 1074580.58 921069.07 1074580.58 46891593.84
54 STN ESS 0.00 1.00000 0.00 STN_ESS 0.00 0.00 0.00 .
55 STN LAM 4273955.04 0.76100 3018.85 STN_LAM 4701350.54 0.00 4701350.54 360935113.7
56 STS ESS 0.00 1.00000 0.00 STS_ESS 0.00 0.00 0.00 .
57 STS LAM 555818.74 0.46234 50901.83 STS_LAM 583609.68 500236.87 583609.68 23586050.24
58 VDO ESS 2206729.31 0.80000 269.71 VDO_ESS 2537738.71 2206729.31 2537738.71 35846491.57
59 VDO LAM 2759217.12 0.76012 81.32 VDO_LAM 3173099.69 2759217.12 3173099.69 7113833.91
60 VDO UNU 0.00 1.00000 0.00 VDO_UNU 0.00 0.00 0.00 .
Spends_change
1 30%
2 -30%
3 30%
4 -30%
5 .
6 .
7 .
8 .
9 10%
10 10%
11 -44%
12 10%
13 .
14 10%
15 10%
16 -5%
17 60%
18 5%
19 2%
20 5%
21 5%
22 -30%
23 -30%
24 -30%
25 -20%
26 20%
27 20%
28 20%
29 -20%
30 -20%
31 20%
32 -10%
33 .
34 5%
35 -10%
36 10%
37 10%
38 10%
39 10%
40 10%
41 10%
42 10%
43 .
44 10%
45 15%
46 0.00%
47 15%
48 -20%
49 -20%
50 -20%
51 .
52 10%
53 5%
54 .
55 10%
56 .
57 5%
58 15%
59 15%
60 .
Upvotes: 1
Views: 4887
Answers (1)
LyzandeR
Reputation: 37879
Took me a while to debug but I made it. First of all I made a few improvements to your functions to make them run faster:
library(Rsolnp)
#first thing to do is calculate the sum of data bounded before subsetting
#the Data_bounded data.frame below
total_spend<-sum(Data_bounded$spend)
#those rows with alpha=0 will always have spend=0 since you maximize.
#So, you don't need them in the calculation. You only cause it
#to be more slow and increase the number of parameters
#when you know that these will always have spend=0
#Therefore I am excluding them
Data_bounded <- na.omit(Data_bounded[Data_bounded$alpha > 0,])
IMPORTANT
Another reason you need to exclude the rows where alpha == 0 is that in many cases you get the lower bound equal to zero and also the upper bound equal to zero i.e. the spend variable cannot take any number and produces NANs and cause the algorithm to fail (as an example see the last row of your data.frame).
alpha<-Data_bounded$alpha
beta<-Data_bounded$beta
#this is what causes big trouble (see below in the answer)
LB_data<-Data_bounded$Lower
UB_data<-Data_bounded$Upper
#I simplified your objective function. Using a for-loop is
#100 times slower
fn1<-function(spend)
{
-sum(alpha * spend^beta)
}
#the constraint equation
eq1<-function(spend)
{
sum(spend)
}
The above settings will cause a non-convergence simply because the range of the lower and upper bounds for each of the parameters are dramatically different:
> summary(LB_data)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0 105400 563800 1045000 1433000 7674000
> summary(UB_data)
Min. 1st Qu. Median Mean 3rd Qu. Max.
81490 529000 1089000 3367000 3206000 37910000
Using these makes the hessian estimation impossible to calculate and the algorithm does not converge:
> solnp(rep(0.5,40),fun=fn1,LB=LB_data,UB=UB_data,eqfun=eq1,eqB=total_spend)
solnp-->The linearized problem has no feasible
solnp-->solution. The problem may not be feasible.
Iter: 1 fn: -504199815.6409 Pars: 4667300.58566 1428556.94211 1769551.34256 0.50002 0.50002 0.50002 0.50002 0.50002 0.50002 920038.05845 1423762.10411 8943897.27854 506400.51947 655335.69677 1427283.02003 1424483.55152 892622.98201 142709.47599 984386.60168 3367895.50511 472422.09101 107401.43750 1010577.85985 57016.73234 60975.13135 1093896.28870 292081.76334 206116.88396 999731.21048 3128070.24034 115116.37073 1248989.03752 111401.37452 5138141.05109 0.50002 1065101.70687 0.50002 508212.67696 2209897.98640 3164060.56280
solnp--> Solution not reliable....Problem Inverting Hessian.
The only solution is to provide good boundaries for your parameters (the sum of estimated spend will still obviously be equal to total_spend so this is for sure correct). And this will work:
#so just set all the boundaries between 0 and Inf.
#solnp will easily be able to find the optimal parameters even though
#the search space is big.
#this took a few seconds on my laptop
opti<-solnp(rep(0.5,40),fun=fn1,LB=rep(0,40),UB=rep(Inf,40),eqfun=eq1,eqB=total_spend)
And you got your convergence and optimal values!!
> opti
$pars
[1] 5.214072e+05 5.457175e+05 1.226056e+02 2.459264e+04 6.727744e+06 4.592510e+02 1.144198e+07 3.482826e+05 2.222721e+07 1.146275e+06 7.594007e+06
[12] 1.338454e+05 1.903417e+05 3.010091e+01 3.705320e+06 1.092090e+00 4.652869e+05 2.935163e+01 3.307376e+03 2.073538e+02 4.628795e+04 1.267862e+03
[23] 4.366613e+02 1.796496e+06 1.396554e+03 3.488297e+06 7.859770e+05 2.858621e+06 1.218379e+02 2.200433e+05 3.241032e+02 3.490860e+04 1.837521e+00
[34] 3.846448e+04 1.237474e+07 1.134505e+07 2.552423e+07 5.201814e+06 1.525503e+06 3.634506e+00
$convergence
[1] 0
EDIT
Now that you showed me the SAS output I can see that you can have a spend amount greater than zero even if the alpha is zero. Keeping this in mind and only excluding cases where the upper and lower limits are zero (spend for these is implied to be zero - same as in the SAS output) this works and converges and uses the actual upper and lower boundaries:
library(Rsolnp)
#first thing to do is calculate the sum of data bounded before subsetting
#the Data_bounded data.frame below
total_spend<-sum(Data_bounded$spend)
Data_bounded <- na.omit(Data_bounded[!(Data_bounded$Lower == 0 & Data_bounded$Upper == 0), ])
alpha<-Data_bounded$alpha
beta<-Data_bounded$beta
#this is what causes big trouble (see below in the answer)
LB_data<-Data_bounded$Lower
UB_data<-Data_bounded$Upper
#I simplified your objective function. Using a for-loop is
#100 times slower
fn1<-function(spend)
{
-sum(alpha * spend^beta)
}
#the constraint equation
eq1<-function(spend)
{
sum(spend)
}
opti<-solnp(rep(0.5,40),fun=fn1,LB=rep(0,40),UB=rep(Inf,40),eqfun=eq1,eqB=total_spend)
Output:
> opti
$pars
[1] 1.863565e+06 8.537037e-02 4.980553e+06 5.078493e+04 1.756843e+07 2.178128e+07 2.455468e+07 1.681695e+07 7.889301e+06 2.746788e+07
[11] 8.545552e-02 5.534938e+05 7.639648e+06 6.134385e+05 7.901136e+05 2.876899e+06 8.545261e-02 8.545261e-02 2.943462e+03 1.178904e+02
[21] 9.985594e+05 9.349150e+02 5.105370e+05 8.303612e+02 8.545261e-02 3.959562e+05 8.545261e-02 8.545261e-02 9.386337e+05 8.545257e-02
[31] 1.607354e+02 2.848127e+04 2.931769e+02 3.296281e+04 8.193187e+05 3.486277e+06 2.300425e+03 1.352874e+05 2.158217e+04 8.545291e-02
$convergence
[1] 0
Upvotes: 3
Related Questions
- R solver optimization
- Optimization R code - Rcpp
- Problems with convergency when using solnp function in R
- Rsolnp: In cbind(temp, funv) : number of rows of result is not a multiple of vector length (arg 1)
- Warning in Rsolnp : Solution not reliable because unable to Invert Hessian
- Optimization using package "nloptr"
- R optimization with optim
- nonlinear optimizaion in R
- Optimisation in R
- Optimisation issue in R