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
-------------- next part --------------
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")))
-------------- next part --------------
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")))
-------------- next part --------------
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"))
-------------- next part --------------
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)
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
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