[R] Fwd: Matrix Constraints in R Optim
Jeff Newmiller
jdnewmil at dcn.davis.ca.us
Tue Jun 21 19:18:50 CEST 2016
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