[R] Fwd: Matrix Constraints in R Optim
Priyank Dwivedi
dpriyank23 at gmail.com
Tue Jun 21 22:13:35 CEST 2016
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
More information about the R-help
mailing list