Skip to content

Matrix Constraints in R Optim

6 messages · Priyank Dwivedi, Jeff Newmiller

#
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,
PD
#
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.
2 days later
#
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:

  
    
1 day later
#
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:

            
---------------------------------------------------------------------------
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:

  
    
#
Not my problem. You are the one applying constraints.