Matrix Constraints in R Optim
Thank you Jeff. It seems to definitely solve it but the "total_error" is very high. Around 399. I also tried with the method = "L-BFGS-B". Still the error is around 399. How can we reduce it? Priyank
On 21 June 2016 at 12:18, Jeff Newmiller <jdnewmil at dcn.davis.ca.us> wrote:
The size of this request is a bit big for this list.
I think you need the constrOptim function to achieve this constraint. See
reproducible example below (no contributed packages needed):
#-----
my.data.matrix.inj <- structure(c(284.6624, 284.6743, 284.6771, 284.6746,
284.6664, 284.6516, 284.6283, 284.5931, 555.1354, 555.0648, 555.0361,
2717.121, 2716.909, 2716.857, 3537.007, 3537.209, 454.2328, 454.2205,
454.2086, 1297.769, 1297.827, 1386.995, 2040.08, 2040.237, 1074.394,
1409.096, 1187.767, 1453.882, 1149.305, 1329.487, 1376.219, 1881.046,
1538.514, 1002.312, 612.8742, 1373.664, 1424.084, 1352.598, 1479.259,
767.9471, 1277.077, 1477.096, 1383.378, 1398.408, 1353.671, 882.6216,
1399.007, 1159.061, 1507.469, 1089.506, 1642.942, 1799.764, 1873.927,
2145.548, 2017.962, 1993.64, 2221.32, 2123.962, 2463.256, 2405.041,
2404.414, 2438.734, 2638.787, 2616.91, 2346.845, 2852.143, 2942.838,
3140.032, 762.2396, 1720.488, 1789.752, 371.4107, 1225.91, 1686.064,
1652.747, 1724.248, 1655.486, 1552.557, 1870.383, 1807.614, 1498.599,
1376.45, 1453.844, 1441.684, 1363.064, 1066.156, 1365.101, 1358.903,
1288.348, 610.3185, 532.7502, 1573.272, 1768.713, 1781.086, 1747.261,
1977.336, 1904.75, 1538.454, 1678.361, 1774.035, 1495.381, 1285.172,
1511.251, 1627.114, 1626.432, 1579.333, 1574.744, 1435.232, 2135.695,
2031.769, 2350.99, 2562.418, 2515.922, 2709.281, 1824.588, 1824.665,
1824.682, 1824.666, 1824.613, 1824.519, 1824.37, 1824.144, 1367.973,
1367.799, 1367.728, 626.0895, 626.0406, 626.0286, 299.3024, 299.3194,
1420.26, 1420.222, 1420.185, 1626.06, 1626.133, 1181.016, 1067.529,
1067.611, 1346.783, 1286.029, 1669.494, 1469.061, 1571.632, 1369.969,
1342.855, 1635.875, 1769.014, 1876.71, 1794.846, 1658.31, 1526.607,
1676.101, 1705.561, 1641.514, 1605.627, 1298.534, 1591.755, 1611.691,
1571.183, 1584.321, 1572.948, 1532.965, 1524.934, 1534.853, 1538.834,
1463.963, 1462.23, 1420.739, 1447.045, 1406.715, 1419.408, 1478.69,
1273.244, 1262.34, 1165.642, 787.8699, 657.2443, 617.5942, 672.4419,
562.5458, 600.0635, 553.3339, 581.2515, 686.7953, 448.5355, 1967.524,
968.7045, 1253.422, 1417.029, 1348.352, 607.6661, 795.2877, 1122.037,
951.7014, 1218.465, 1452.847, 1708.894, 1789.318, 1774.066, 1730.023,
1792.384, 1647.639, 1532.214, 1398.604, 1456.599, 1405.635, 1341.6,
1384.088, 1547.139, 1480.687, 1527.453, 1541.885, 1348.729, 1359.007,
1093.668, 1078.121, 1202.416, 895.9857, 1175.532, 1010.464, 967.2054,
851.1081, 740.4431, 930.6541, 1057.503, 1036.018, 1250.418, 1382.047,
278.5883, 278.6001, 278.6027, 278.6003, 278.5922, 278.5778, 278.555,
278.5205, 922.1713, 922.054, 922.0063, 967.21, 967.1343, 967.1157, 774.6002,
774.6443, 772.0591, 772.0382, 772.018, 870.8308, 870.8698, 904.8117,
825.425, 792.9547, 876.54, 882.7752, 681.8339, 775.945, 1081.869, 928.0758,
921.4498, 1079.74, 795.1276, 810.2282, 835.9764, 825.9167, 825.2587,
943.9789, 745.8108, 709.2183, 718.3409, 656.4478, 553.8104, 682.1406,
863.1352, 837.0597, 850.8278, 789.4566, 827.334, 813.5239, 723.0217,
808.3031, 871.0251, 1023.663, 1008.41, 1118.704, 1113.178, 907.1134,
726.4997, 1064.354, 1208.275, 1269.964, 1226.312, 834.8596, 952.5037,
1019.817, 922.9584, 886.3052, 898.9753, 868.3756, 869.4521, 105.3649,
407.1053, 136.8827, 722.5133, 841.0006, 706.9567, 542.9826, 198.147,
233.6965, 114.3593, 252.4854, 284.9101, 418.044, 215.6109, 543.6895,
654.181, 927.2443, 896.0264, 822.9401, 878.3534, 692.4314, 738.8477,
984.3605, 1069.655, 1022.925, 1002.807, 850.6902, 991.8134, 1034.01,
1148.745, 1142.539, 1163.838, 1275.52, 1145.691, 1460.11, 1377.891,
1306.395, 1304.617, 1278.456, 1378.95, 1374.073, 1449.972, 1184.909,
270.0509, 270.0623, 270.0648, 270.0624, 270.0547, 270.0407, 270.0186,
269.9851, 631.4337, 631.3534, 631.3207, 607.871, 607.8235, 607.8118,
570.1067, 570.1392, 471.9973, 471.9845, 471.9722, 374.8601, 374.8769,
482.6559, 509.4759, 259.6612, 601.047, 612.4909, 599.3603, 368.4525,
541.0823, 637.376, 572.6561, 520.8604, 602.978, 508.6731, 518.9494,
559.4774, 583.3226, 665.8262, 675.3377, 604.7722, 619.4575, 567.0582,
700.1987, 680.9487, 720.6385, 697.012, 662.4166, 683.2136, 659.8345,
667.4672, 707.6854, 743.7268, 858.9992, 832.3246, 779.6216, 698.0973,
703.4314, 791.7886, 726.9083, 854.6981, 834.7772, 832.3445, 812.7689,
727.6645, 652.1965, 826.9865, 849.4389, 811.6799, 850.7483, 832.3735,
819.6655, 1042.436, 720.7501, 952.0648, 1195, 848.0734, 976.9899, 1112.395,
1113.345, 1153.728, 805.5801, 646.0727, 617.1312, 791.8318, 847.233,
683.816, 724.7269, 911.1725, 827.3728, 995.0048, 800.6775, 879.0817,
972.6709, 799.3595, 1029.595, 1007.769, 852.9899, 837.8101, 941.9149,
982.4396, 979.9702, 967.2394, 937.1133, 960.9035, 908.2497, 996.8404,
1190.648, 1202.747, 1350.496, 1267.897, 1132.526, 1055.183, 799.7894,
639.9702, 769.6429, 769.6754, 769.6827, 769.676, 769.6537, 769.6139,
769.551, 769.4556, 499.9228, 499.8593, 499.8334, 1051.619, 1051.537,
1051.517, 1017.837, 1017.895, 787.5231, 787.5018, 787.4812, 127.3492,
127.3549, 240.9772, 248.1084, 400.2578, 663.3332, 986.2067, 936.059,
1061.159, 849.0998, 884.3383, 1183.185, 1208.31, 981.9471, 1076.72,
1124.325, 1008.958, 780.2723, 692.6738, 1044.181, 804.3527, 664.2988,
713.3538, 768.6463, 791.4983, 1408.636, 1460.505, 1331.472, 1436.979,
1223.143, 1192.528, 1165.123, 1187.325, 889.4554, 1755.404, 1539.565,
1367.623, 1197.647, 1204.832, 1253.376, 1064.125, 1221.669, 1063.684,
1029.96, 941.9225, 953.305, 1135.038, 995.6816, 1202.049, 1179.09, 1238.77,
1252.872, 195.4976, 796.9503, 1409.675, 2215.336, 1971.793, 1372.014,
1194.094, 990.832, 1240.13, 1272.831, 1110.265, 1083.954, 1277.695,
1224.066, 1216.931, 1036.133, 1275.89, 650.2736, 493.1569, 443.461,
457.3099, 492.6304, 514.841, 490.7231, 505.4785, 567.1318, 544.3971,
547.5244, 528.4097, 662.0999, 964.6831, 1006.148, 1102.357, 1207.62,
1272.277, 1173.155, 1125.227, 1039.502, 1074.456, 1146.245, 1429.14,
1246.974, 1215.329)
, .Dim = c(114L, 5L)
, .Dimnames = list(NULL, c("I1", "I2", "I3", "I4", "I5")))
my.data.matrix.prod <- structure(c(2916.28, 1893.82, 1446.496, 1223.643,
1093.515, 1027.691, 1025.575, 1069.484, 1350.653, 1383.106, 1404.12,
3229.087, 3287.819, 3292.214, 3949.526, 3934.924, 1344.882, 1276.475,
1281.724, 2080.675, 2170.162, 2204.06, 2733.114, 2709.72, 1906.547,
2226.197, 2147.538, 2396.16, 2170.339, 2295.214, 2325.382, 2863.881,
2633.29, 2191.615, 1823.576, 2462.448, 2472.716, 2426.248, 2558.359,
1898.222, 2311.003, 2405.334, 2359.773, 2406.227, 2404.66, 2005.993,
2470.426, 2262.771, 2564.288, 2187.93, 2672.702, 2817.843, 2886.186,
3159.216, 3071.983, 3038.874, 3232.614, 3153.618, 3396.065, 3337.943,
3314.298, 3228.766, 3312.479, 3214.223, 2943.438, 3374.134, 3471.613,
3649.256, 1494.396, 2318.848, 2353.137, 1392.929, 2017.725, 2497.875,
2650.34, 2772.884, 2503.756, 2341.685, 2665.939, 2603.909, 2361.046,
2307.904, 2466.254, 2545.271, 2505.55, 2239.917, 2518.568, 2521.566,
2398.009, 1700.699, 1570.964, 2475.785, 2666.551, 2696.887, 2733.822,
2956.056, 2906.461, 2566.767, 2639.433, 2717.689, 2399.816, 2175.098,
2405.237, 2461.575, 2513.077, 2476.729, 2467.291, 2303.615, 2898.341,
2858.363, 3200.795, 3426.61, 3443.722, 3647.533, 195.3348, 176.5879,
161.8616, 147.6775, 132.3667, 116.3203, 100.9762, 90.91395, 102.5056,
111.2312, 119.294, 139.5639, 148.0501, 154.4379, 162.0608, 166.5477,
150.7256, 143.1064, 137.4059, 131.9734, 127.8249, 129.6863, 136.1022,
121.6995, 131.2575, 144.92, 150.0162, 140.1022, 146.21, 156.451, 158.8145,
162.6809, 164.6031, 156.6059, 150.636, 155.6411, 158.7302, 166.222,
171.0211, 162.1327, 161.2135, 156.3216, 162.0996, 166.6428, 175.9184,
176.8375, 178.7133, 178.7524, 179.3178, 176.1973, 180.4867, 187.3193,
199.2127, 209.983, 210.1795, 203.8254, 201.1218, 203.3554, 199.8475,
209.8946, 215.1455, 215.6018, 213.2702, 199.2345, 185.3278, 197.2057,
205.0727, 207.3002, 193.6611, 194.2139, 193.8643, 193.2228, 177.8776,
191.4582, 231.8191, 227.0726, 224.6594, 229.7895, 230.8227, 234.7284,
206.1662, 179.8467, 167.3609, 179.5722, 188.3897, 180.9705, 182.7036,
202.3105, 200.8232, 203.9204, 189.2181, 192.9931, 204.6493, 199.082,
215.5948, 223.7031, 213.8644, 202.6964, 208.5682, 216.1876, 217.9815,
217.007, 217.463, 221.4278, 218.8876, 228.6546, 247.8913, 255.3423,
274.8202, 276.3341, 269.6512, 262.6747, 239.2566, 213.2598, 196.0692,
179.3542, 174.4489, 179.1992, 193.516, 219.7416, 261.9235, 307.7595,
339.0413, 349.1725, 355.6877, 355.0119, 353.4153, 351.7466, 334.9937,
315.7924, 338.9163, 353.4399, 367.3095, 370.7577, 368.3222, 338.1546,
309.5753, 302.9909, 343.3383, 390.3582, 442.8708, 467.6517, 475.7294,
463.2386, 475.4719, 512.9818, 525.4725, 546.562, 555.4177, 539.306,
499.4974, 483.7216, 504.7977, 493.012, 470.5119, 433.9357, 442.8588,
456.0057, 512.4643, 550.1924, 558.0298, 564.4106, 550.0839, 538.8026,
530.6313, 523.3772, 495.919, 552.271, 570.1813, 559.772, 539.137, 531.285,
511.5488, 489.0468, 483.7139, 434.6737, 391.9633, 353.6852, 341.9161,
345.1014, 337.9316, 347.3225, 351.3463, 368.2297, 356.9464, 385.1567,
367.8657, 433.608, 567.5147, 609.7797, 502.3615, 441.0474, 421.461,
421.8376, 446.8344, 468.2683, 499.6648, 542.9484, 556.0471, 560.2142,
552.7231, 561.0404, 498.0082, 435.9498, 406.5746, 388.8749, 379.8109,
384.6039, 402.4961, 406.7456, 417.0511, 418.7817, 404.1004, 396.1866,
381.1434, 398.5426, 424.4879, 419.1766, 448.4539, 459.9056, 450.9682,
429.8293, 402.7214, 409.8873, 434.7366, 470.5877, 491.6042, 505.3956,
2379.811, 1683.061, 1348.136, 1183.511, 1096.342, 1063.209, 1083.307,
1137.872, 1698.039, 1777.531, 1824.798, 1990.391, 2049.531, 2094.436,
1982.723, 1974.184, 1931.659, 1916.844, 1909.946, 1859.683, 1768.624,
1733.896, 1644.874, 1566.683, 1802.985, 2026.399, 2002.01, 2095.246,
2341.096, 2261.631, 2337.393, 2534.549, 2322.27, 2333.124, 2367.872,
2336.886, 2235.952, 2284.032, 2240.623, 2144.319, 2069.301, 1946.398,
1910.047, 2043.538, 2433.989, 2556.197, 2578.596, 2568.367, 2534.411,
2478.992, 2395.445, 2468.419, 2459.169, 2850.418, 2889.555, 2898.655,
2802.34, 2630.252, 2451.473, 2667.949, 2813.618, 2777.629, 2657.484,
2226.911, 2225.193, 2366.028, 2296.084, 2321.493, 2335.952, 2351.763,
2347.124, 1576.808, 1743.489, 1847.67, 2869.197, 3040.427, 2686.383,
2409.108, 2028.894, 2091.013, 1932.818, 1923.021, 1920.55, 2171.707,
2086.544, 2304.832, 2344.914, 2685.513, 2448.996, 2252.836, 2147.083,
1971.758, 2033.175, 2196.932, 2353.921, 2357.346, 2326.293, 2178.859,
2293.083, 2341.083, 2452.64, 2557.318, 2645.425, 2778.83, 2744.436,
3066.146, 3070.198, 3004.751, 3008.488, 2991.268, 3074.413, 3159.114,
3126.801, 2823.369)
, .Dim = c(114L, 4L)
, .Dimnames = list(NULL, c("Q1", "Q2", "Q3", "Q4")))
my.data.matrix.time <- structure(c(1, 1.944202, 3.803123, 6.203458,
9.420446, 14.03878, 21.35927, 30.4375, 44.67685, 52.77593, 60.875, 76.09375,
83.70312, 91.3125, 104.9416, 121.75, 136.9688, 144.5781, 152.1875, 167.4062,
182.625, 213.0625, 243.5, 273.9375, 304.375, 334.8125, 365.25, 395.6875,
426.125, 456.5625, 487, 517.4375, 547.875, 578.3125, 608.75, 639.1875,
669.625, 700.0625, 730.5, 760.9375, 791.375, 821.8125, 852.25, 882.6875,
913.125, 943.5625, 974, 1004.438, 1034.875, 1065.312, 1095.75, 1126.188,
1156.625, 1187.062, 1217.5, 1247.938, 1278.375, 1308.812, 1339.25, 1369.688,
1400.125, 1430.562, 1461, 1491.438, 1521.875, 1552.312, 1582.75, 1613.188,
1643.625, 1674.062, 1704.5, 1734.938, 1765.375, 1795.812, 1826.25, 1856.688,
1887.125, 1917.562, 1948, 1978.438, 2008.875, 2039.312, 2069.75, 2100.188,
2130.625, 2161.062, 2191.5, 2221.938, 2252.375, 2282.812, 2313.25, 2343.688,
2374.125, 2404.562, 2435, 2465.438, 2495.875, 2526.312, 2556.75, 2587.188,
2617.625, 2648.062, 2678.5, 2708.938, 2739.375, 2769.812, 2800.25, 2830.688,
2861.125, 2891.562, 2922, 2952.438, 2982.875, 3013.312)
, .Dim = c(114L, 1L)
, .Dimnames = list( NULL, "time"))
my.data.var <- c( 10, 0.25, 0.25, 0.25, 0.25, 0.25
, 10, 0.25, 0.25, 0.25, 0.25, 0.25
, 10, 0.25, 0.25, 0.25, 0.25, 0.25
, 10, 0.25, 0.25, 0.25, 0.25, 0.25
)
my.data.qo <- c( 5990, 150, 199, 996 ) #Pre-Waterflood Production
my.data.timet0 <- 0 # starting condition for time
#FUNCTION
Qjk.Cal.func <- function( my.data.timet0
, my.data.qo
, my.data.matrix.time
, my.data.matrix.inj
, my.data.matrix.prod
, my.data.var
, my.data.var.mat
)
{
qjk.cal.matrix <- matrix(
, nrow = nrow( my.data.matrix.prod )
, ncol = ncol( my.data.matrix.prod )
)
count <- 1
number <- 1
# loop through all PROD wells columns
for ( colnum in 1:ncol( my.data.matrix.prod ) ) {
# "sum" is a very bad choice of variable name
# as it is a commonly-used base function
sum <- 0 # this initialization is redundant see below
#loop through all the rows
for( row in 1:nrow( my.data.matrix.prod ) ) {
sum <- 0 # most frequent re-initialization
deltaT <- 0
expo <- 0
#loop through all the injector columns to get the PRODUCT SUM
for( column in 1:ncol( my.data.matrix.inj ) ) {
sum <- ( sum
+ my.data.matrix.inj[ row, column ]
* my.data.var.mat[ colnum, number+column ]
)
}
if ( count < 2 ) {
deltaT <- my.data.matrix.time[ row ]
} else {
deltaT <- ( my.data.matrix.time[ row ]
- my.data.matrix.time[ row - 1 ]
)
}
expo <- exp( -deltaT / my.data.var.mat[ colnum, 1 ] )
# change here too
if ( count < 2 ) {
qjk.cal.matrix[ row, colnum ] <-
my.data.qo[ colnum ] * expo + ( 1 - expo ) * sum
} else {
qjk.cal.matrix[ row, colnum ] <-
( qjk.cal.matrix[ row - 1, colnum ] * expo
+ ( 1 - expo ) * sum
)
}
count <- count + 1
}
count <- 1
}
# RETURN CALCULATED MATRIX TO THE ERROR FUNCTION
qjk.cal.matrix
}
Error.func <- function( my.data.var ) {
#First convert vector(my.data.var) to MATRIX
# and send it to calculate new MATRIX
my.data.var.mat <- matrix( my.data.var
, nrow = ncol( my.data.matrix.prod )
, ncol = ncol( my.data.matrix.inj ) + 1
, byrow = TRUE
)
Calc.Qjk.Value <- Qjk.Cal.func( my.data.timet0
, my.data.qo
, my.data.matrix.time
, my.data.matrix.inj
, my.data.matrix.prod
, my.data.var
, my.data.var.mat
)
#FIND DIFFERENCE BETWEEN CAL. MATRIX AND ORIGINAL MATRIX
diff.values <- my.data.matrix.prod - Calc.Qjk.Value
#sum of square root of the diff
Error <- ( ( colSums( ( diff.values^2 )
, na.rm = FALSE
, dims = 1
)
) / nrow( my.data.matrix.inj )
)^0.5
#print(paste(Error))
# total avg error
Error_total <- sum( Error
, na.rm=FALSE
) / ncol( my.data.matrix.prod )
Error_total
}
n <- ncol( my.data.matrix.prod )
m <- ncol( my.data.matrix.inj )
k <- ( 1 + n ) * m
ciA <- numeric( k )
uiA <- array( 0, dim = c( m+1, n, k ) )
ia <- 0
#Aoff2 <- (m+1) * n
for ( i in seq.int( m ) ) {
ia <- ia + 1L
# sum of columns <= 1
uiA[ i+1, , i ] <- -1
ciA[ ia ] <- -1
}
for ( i in ( 1 + seq.int( m ) ) ) {
for ( j in seq.int( n ) ) {
ia <- ia + 1L
# elements > 0
uiA[ i, j, ia ] <- 1
ciA[ ia ] <- 0
}
}
uiA <- matrix( uiA, nrow = k, byrow = TRUE )
my.data.varA <- c( 10, 0.25, 0.25, 0.25, 0.25, 0.25
, 10, 0.25, 0.25, 0.25, 0.25, 0.25
, 10, 0.25, 0.25, 0.25, 0.25, 0.25
, 10, 0.24, 0.24, 0.24, 0.24, 0.24
)
# interior?
all( uiA %*% my.data.var1 - (ciA) > 0 )
sols <- constrOptim( my.data.varA
, Error.func
, NULL
, ui = uiA
, ci = ciA
, method = "SANN"
)
# meets constraint?
all( uiA %*% sols$par - (ciA) >= 0 )
sols
##########################
On Sun, 19 Jun 2016, Priyank Dwivedi wrote:
All,
Here are the dput files of the input data to the code.
Thanks for any advice.
I am adding the entire code below too just in case.
file <- file.path("Learning R","CRM_R_Ver4.xlsx")
file
my.data <- readWorksheetFromFile(file,sheet=1,startRow=1)
str(my.data) # DATA FRAME
my.data.matrix.inj <- as.matrix(my.data) #convert DATA FRAME to MATRIX
my.data.matrix.inj
dput(my.data.matrix.inj,"my.data.matrix.inj.txt")
my.data.2 <- readWorksheetFromFile(file,sheet=2,startRow=1)
str(my.data.2) # DATA FRAME
my.data.matrix.time <- as.matrix(my.data.2) #convert DATA FRAME to MATRIX
my.data.matrix.time
dput(my.data.matrix.time,"my.data.matrix.time.txt")
my.data <- readWorksheetFromFile(file,sheet=3,startRow=1)
str(my.data) # DATA FRAME
my.data.matrix.prod <- as.matrix(my.data) #convert DATA FRAME to MATRIX
my.data.matrix.prod
dput(my.data.matrix.prod,"my.data.matrix.prod.txt")
# my.data.var <- vector("numeric",length = 24)
# my.data.var
my.data.var <- c(10,0.25,0.25,0.25,0.25,0.25,
10,0.25,0.25,0.25,0.25,0.25,
10,0.25,0.25,0.25,0.25,0.25,
10,0.25,0.25,0.25,0.25,0.25)
my.data.var
dput(my.data.var,"my.data.var.txt")
my.data.qo <- c(5990,150,199,996) #Pre-Waterflood Production
my.data.timet0 <- 0 # starting condition for time
#FUNCTION
Qjk.Cal.func <- function(my.data.timet0,my.data.qo,my.data.matrix.time,
my.data.matrix.inj,
my.data.matrix.prod,my.data.var,my.data.var.mat)
{
qjk.cal.matrix <- matrix(,nrow = nrow(my.data.matrix.prod),
ncol=ncol(my.data.matrix.prod))
count <- 1
number <- 1
for(colnum in 1:ncol(my.data.matrix.prod)) # loop through all PROD
wells columns
{
sum <-0
for(row in 1:nrow(my.data.matrix.prod)) #loop through all the rows
{
sum <-0
deltaT <-0
expo <-0
for(column in 1:ncol(my.data.matrix.inj)) #loop through all
the injector columns to get the PRODUCT SUM
{
sum = sum +
my.data.matrix.inj[row,column]*my.data.var.mat[colnum,number+column]
}
if(count<2)
{
deltaT<- my.data.matrix.time[row]
}
else
{deltaT <- my.data.matrix.time[row]-my.data.matrix.time[row-1]}
expo <- exp(-deltaT/my.data.var.mat[colnum,1])
# change here too
if(count<2)
{
qjk.cal.matrix[row,colnum] = my.data.qo[colnum]*expo + (1-expo)*sum
}
else
{
qjk.cal.matrix[row,colnum]=qjk.cal.matrix[row-1,colnum]*expo +
(1-expo)*sum
}
count <- count+1
}
count <-1
}
qjk.cal.matrix # RETURN CALCULATED MATRIX TO THE ERROR FUNCTION
}
# ERROR FUNCTION - FINDS DIFFERENCE BETWEEN CAL. MATRIX AND ORIGINAL
MATRIX. Miminize the Error by changing my.data.var
Error.func <- function(my.data.var)
{
#First convert vector(my.data.var) to MATRIX aand send it to
calculate new MATRIX
my.data.var.mat <- matrix(my.data.var,nrow =
ncol(my.data.matrix.prod),ncol = ncol(my.data.matrix.inj)+1,byrow =
TRUE)
Calc.Qjk.Value <-
Qjk.Cal.func(my.data.timet0,my.data.qo,my.data.matrix.time,
my.data.matrix.inj,
my.data.matrix.prod,my.data.var,my.data.var.mat)
diff.values <- my.data.matrix.prod-Calc.Qjk.Value #FIND
DIFFERENCE BETWEEN CAL. MATRIX AND ORIGINAL MATRIX
Error <- ((colSums ((diff.values^2), na.rm = FALSE, dims =
1))/nrow(my.data.matrix.inj))^0.5 #sum of square root of the diff
print(paste(Error))
Error_total <- sum(Error,na.rm=FALSE)/ncol(my.data.matrix.prod) #
total avg error
Error_total
}
# OPTIMIZE
sols<-optim(my.data.var,Error.func,method="L-BFGS-B",upper=c(Inf,1,1,1,1,1,Inf,1,1,1,1,1,Inf,1,1,1,1,1,Inf,1,1,1,1,1),
lower=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))
sols
On 17 June 2016 at 16:55, Jeff Newmiller <jdnewmil at dcn.davis.ca.us> wrote:
Your code is corrupt because you failed to send your email in plain text format. You also don't appear to have all data needed to reproduce the problem. Use the dput function to generate R code form of a sample of your data. -- Sent from my phone. Please excuse my brevity. On June 17, 2016 1:07:21 PM PDT, Priyank Dwivedi <dpriyank23 at gmail.com> wrote:
By mistake, I sent it earlier to the wrong address.
---------- Forwarded message ----------
From: Priyank Dwivedi <dpriyank23 at gmail.com>
Date: 17 June 2016 at 14:50
Subject: Matrix Constraints in R Optim
To: r-help-owner at r-project.org
Hi,
Below is the code snippet I wrote in R:
The basic idea is to minimize error by optimizing set of values (in this
scenario 12) in the form of a matrix. I defined the matrix elements as
vector "*my.data.var" * and then stacked it into a matrix called
"*my.data.var.mat"
in the error function. *
The only part that I can't figure out is "what if the column sum in
the *my.data.var.mat
needs to be <=1"; that's the constraint/s.. Where do I introduce it in
the
OPTIM solver or elsewhere?*
*my.data.matrix.inj* <- as.matrix(my.data) #convert DATA FRAME to
MATRIX
my.data.matrix.inj
*my.data.matrix.time* <- as.matrix(my.data.2) #convert DATA FRAME to
MATRIX
my.data.matrix.time
*my.data.matrix.prod* <- as.matrix(my.data) #convert DATA FRAME to
MATRIX
my.data.matrix.prod
*my.data.var* <-
c(2,0.8,0.5,0.2,0.2,0.1,10,0.01,0.02,0.2,0.1,0.01,2,0.8,0.5,0.2,0.2,0.1,10,0.01,0.02,0.2,0.1,0.01,2,0.8,0.5,0.2,0.2,0.1,10,0.01,0.02,0.2,0.1,0.01)
my.data.var
*my.data.qo* <- c(5990,150,199,996) #Pre-Waterflood Production
*my.data.timet0* <- 0 # starting condition for time
*#FUNCTIONQjk.Cal.func* <-
function(my.data.timet0,my.data.qo,my.data.matrix.time,
my.data.matrix.inj,
my.data.matrix.prod,my.data.var,my.data.var.mat)
{
qjk.cal.matrix <- matrix(,nrow = nrow(my.data.matrix.prod),
ncol=ncol(my.data.matrix.prod))
count <- 1
number <- 1
for(colnum in 1:ncol(my.data.matrix.prod)) # loop through all PROD
wells columns
{
sum <-0
for(row in 1:nrow(my.data.matrix.prod)) #loop through all the rows
{
sum <-0
deltaT <-0
expo <-0
for(column in 1:ncol(my.data.matrix.inj)) #loop through all the
injector columns to get the PRODUCT SUM
{
sum = sum +
my.data.matrix.inj[row,column]*my.data.var.mat[colnum,number+column]
}
if(count<2)
{
deltaT<- my.data.matrix.time[row]
}
else
{deltaT <- my.data.matrix.time[row]-my.data.matrix.time[row-1]}
expo <- exp(-deltaT/my.data.var.mat[colnum,1]) #
change here too
if(count<2)
{
qjk.cal.matrix[row,colnum] = my.data.qo[colnum]*expo +
(1-expo)*sum
}
else
{
qjk.cal.matrix[row,colnum]=qjk.cal.matrix[row-1,colnum]*expo +
(1-expo)*sum
}
count <- count+1
}
count <-1
}
qjk.cal.matrix # RETURN CALCULATED MATRIX TO THE ERROR FUNCTION
}
*# ERROR FUNCTION* - FINDS DIFFERENCE BETWEEN CAL. MATRIX AND ORIGINAL
MATRIX. Miminize the Error by changing my.data.var
*Error.func* <- function(my.data.var)
{
#First convert vector(my.data.var) to MATRIX aand send it to calculate
new MATRIX
*my.data.var.mat* <- matrix(my.data.var,nrow =
ncol(my.data.matrix.prod),ncol = ncol(my.data.matrix.inj)+1,byrow =
TRUE)
* Calc.Qjk.Value* <-
Qjk.Cal.func(my.data.timet0,my.data.qo,my.data.matrix.time,
my.data.matrix.inj,
my.data.matrix.prod,my.data.var,my.data.var.mat)
diff.values <-
my.data.matrix.prod-Calc.Qjk.Value #FIND DIFFERENCE
BETWEEN CAL. MATRIX AND ORIGINAL MATRIX
Error <- ((colSums ((diff.values^2), na.rm = FALSE, dims =
1))/nrow(my.data.matrix.inj))^0.5 #sum of square root of the diff
print(paste(Error))
Error_total <- sum(Error,na.rm=FALSE)/ncol(my.data.matrix.prod) #
total
avg error
* Error_total*
}
# OPTIMIZE
*optim*(*my.data.var*
,Error.func,method="L-BFGS-B",upper=c(Inf,1,1,1,1,1,Inf,1,1,1,1,1,Inf,1,1,1,1,1,Inf,1,1,1,1,1))
-- Best Regards, Priyank Dwivedi
---------------------------------------------------------------------------
Jeff Newmiller The ..... ..... Go Live...
DCN:<jdnewmil at dcn.davis.ca.us> Basics: ##.#. ##.#. Live Go...
Live: OO#.. Dead: OO#.. Playing
Research Engineer (Solar/Batteries O.O#. #.O#. with
/Software/Embedded Controllers) .OO#. .OO#. rocks...1k
---------------------------------------------------------------------------
Best Regards, Priyank Dwivedi