[R] Fwd: Matrix Constraints in R Optim
Jeff Newmiller
jdnewmil at dcn.davis.ca.us
Tue Jun 21 23:02:55 CEST 2016
Not my problem. You are the one applying constraints.
--
Sent from my phone. Please excuse my brevity.
On June 21, 2016 1:13:35 PM PDT, Priyank Dwivedi <dpriyank23 at gmail.com> wrote:
>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
>>
>---------------------------------------------------------------------------
More information about the R-help
mailing list