Piecewise Cost Representation in Matrix Export

[msdogan]

@jrmerz @qjhart the piecewise cost representation seems to have some bugs in it. We have one extra link in the beginning (when k = 0), and it is constrained. Below, there is a link with piecewise cost HU507 to AGG_COACH. There are two months: July and August. I wonder why it constrains first pieces (k=0) with zero cost. This is true for all links with piecewise representation.

i   j   k   cost    amplitude   lower_bound upper_bound

HU507.2002-07-31    AGG_COACH.2002-07-31    0   0   1   18.49   18.49
HU507.2002-07-31    AGG_COACH.2002-07-31    1   -2450.118243    1   0   1.85
HU507.2002-07-31    AGG_COACH.2002-07-31    2   -2381.098208    1   0   1.84
HU507.2002-07-31    AGG_COACH.2002-07-31    3   -2243.863863    1   0   1.850001
HU507.2002-07-31    AGG_COACH.2002-07-31    4   -1961.753071    1   0   1.849998
HU507.2002-07-31    AGG_COACH.2002-07-31    5   -1641.658741    1   0   1.850001
HU507.2002-07-31    AGG_COACH.2002-07-31    6   -1423.632549    1   0   1.85
HU507.2002-07-31    AGG_COACH.2002-07-31    7   -1208.016126    1   0   1.85
HU507.2002-07-31    AGG_COACH.2002-07-31    8   -979.330264 1   0   1.849999
HU507.2002-07-31    AGG_COACH.2002-07-31    9   -768.2532673    1   0   1.850002
HU507.2002-07-31    AGG_COACH.2002-07-31    10  -337.3858576    1   0   1.84
HU507.2002-07-31    AGG_COACH.2002-07-31    11  0   1   0   0

HU507.2002-08-31    AGG_COACH.2002-08-31    0   0   1   15  15
HU507.2002-08-31    AGG_COACH.2002-08-31    1   -2453   1   0   1.5
HU507.2002-08-31    AGG_COACH.2002-08-31    2   -2369.833984    1   0   1.5
HU507.2002-08-31    AGG_COACH.2002-08-31    3   -2231.628984    1   0   1.51
HU507.2002-08-31    AGG_COACH.2002-08-31    4   -1963.080078    1   0   1.5
HU507.2002-08-31    AGG_COACH.2002-08-31    5   -1643.593099    1   0   1.5
HU507.2002-08-31    AGG_COACH.2002-08-31    6   -1424.586914    1   0   1.5
HU507.2002-08-31    AGG_COACH.2002-08-31    7   -1209.439941    1   0   1.5
HU507.2002-08-31    AGG_COACH.2002-08-31    8   -979.9932453    1   0   1.5
HU507.2002-08-31    AGG_COACH.2002-08-31    9   -769.160034 1   0   1.5
HU507.2002-08-31    AGG_COACH.2002-08-31    10  -335.9533287    1   0   1.5
HU507.2002-08-31    AGG_COACH.2002-08-31    11  0   1   0   0

Aha, I think I got the reason. For HU507 to AGG_COACH for July, our penalty data starts from 18.49. Below this number there is nothing. Because of that, your code constrains it. Below are figures from penalty data DSS.

I think the decision was to extend the beginning and ending penalties. As far as I know, HEC-PRM does that automatically. Am I right @qjhart @mimijenkins1?

[jrmerz]

@msdogan I'm assuming you want to wait on a response to implement this change?

[msdogan]

@jrmerz yes. We discussed this before, but I will bring this up again in the next hobbes meeting.

[msdogan]

PROBLEM:
Extending penalty curve and operating cost curve lines when exporting network matrix, which is done automatically in CALVIN runs by HEC-PRM but for PyVIN we need to do manually.

WHAT HEC-PRM IS DOING:

example: Shasta hydropower storage curve
Original curve is on the right (blue), extended curve is on the left (red)

example: Coachella ag penalty curve
Original curve is on top (blue), extended curve is in bottom (red). It is not shown here but on the x axis it is extended to very large number (10^9)

example: Banks pumping plant operating cost curve
Original curve is on top (blue), extended is in bottom (red). It is also not shown here but on the right hand side, the curve is extended to very large number.

DISCUSSION:

1. I am not sure why hydropower storage curve for Shasta is not extended on the right hand side, but I think it should be extended to x point where y is zero with a slope of ending piece.
2. I don't see a reason why to extend the curve to very large number for Coachella while it is already crossing x axis. The aim here to cross x and y lines.

WHAT WE SHOULD DO:

I don't know how matrix exporting script works but adding a 'if statement' might solve the issue. If curves don't cross axis lines, then extend. And, for operating cost curves, extend it to very large number (max ub) on the right hand side.

Any comments? @jdherman @josue-medellin @mimijenkins1 @qjhart

[mimijenkins1]

Mustafa's suggestions for extending the penalty/cost curves to intersect
the y-axis (left side extension) for the matrix for PyVIN make good sense
and agree with the HEC PRM manual (see Fig 8).

However, for a downward sloping penalty that does not intersect the x-axis
(right side extension, as in the shasta hydropower storage penalty curve),
HEC-PRM manual seems to indicate the extension they do is "a zero cost
slope" (see Fig 8, b & c). If the RED lines you plotted the PD from the
HEC-PRM DSS PD output file, it might be b/c HEC PRM assumes zero cost
slope and this isn't stored in the PD?

Either way, I would tend to agree with extending the negative cost slope
of the last segment on the right to intersect the x-axis if this extension
is needed, for any downward sloping penalty curves.

On Tue, 13 Sep 2016 10:29:50 -0700, Mustafa Dogan
notifications@github.com wrote:

PROBLEM:Extending penalty curve and operating cost curve lines when
exporting network matrix, which is done automatically in CALVIN runs >by
HEC-PRM but for PyVIN we need to do manually.

WHAT HEC-PRM IS DOING:

example: Shasta hydropower storage curve
Original curve is on the right (blue), extended curve is on the left
(red)

example: Coachella ag penalty curve
Original curve is on top (blue), extended curve is in bottom (red). It
is not shown here but on the x axis it is extended to very large >number
(10^9)

example: Banks pumping plant operating cost curve
Original curve is on top (blue), extended is in bottom (red). It is also
not shown here but on the right hand side, the curve is extended to

very large number.

DISCUSSION:
I am not sure why hydropower storage curve for Shasta is not extended on
the right hand side, but I think it should be extended >to x point where
y is zero with a slope of ending piece.
I don't see a reason why to extend the curve to very large number for
Coachella while it is already crossing x axis. The aim here >to cross x
and y lines.

WHAT WE SHOULD DO:

I don't know how matrix exporting script works but adding a 'if
statement' might solve the issue. If curves don't cross axis lines, then

extend. And, for operating cost curves, extend it to very large number
(max ub) on the right hand side.

Any comments? @jdherman @josue-medellin @mimijenkins1 @qjhart
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub, or mute the thread.

Marion W. Jenkins, PhD
Associate Research Engineer
Dept. Civil & Environmental Engineering
Ghausi Hall, Rm 3019
University of California at Davis
Davis, CA 95616
Tel: +1 530 754 6427
Fax: +1 530 752 7872

[jdherman]

Couple of thoughts:

• Are there flags in the data to separate penalty vs. pumping curves? Or would you just have to look at whether the slope is +/-.
• I agree all penalty curves should extend to 0 on the left. This should be easy, code-wise: if the lower bound of the first piecewise piece is not 0, make it 0.
• The right-side extensions are tricker. For example the Shasta one isn't going past 4550 TAF because that's the storage capacity. There are probably lots of similar cases where the values of "U" in the matrix will take precedence over the cost curves.

I think our goal in general is to make sure the full operating range of these reservoirs/pumps/canals/whatever all have costs assigned to them. Anything outside of that range shouldn't matter.

[msdogan]

@mimijenkins1 all red lines are coming from HEC-PRM output DSS file. For hydropower penalty, our case is Fig 8 (c) is HEC-PRM Manual. Looks like the curve is extrapolated to y-axis but not x-axis. Below is paired data (from output DSS) of the same plot above.

@jdherman we don't have any flags for penalty data. As you wrote, it is based on +/- whether it is operating cost or demand penalty. For your third point, yes there is a storage capacity and it makes it a little complicating. Maybe that's why HEC-PRM didn't extend it. Then, if there is a capacity (upper bound), we should extend until upper bound?

I think current code cuts the cost curve at the physical upper bound (capacity) value if it extends more than capacity when creating piecewise pieces. Am I right @jrmerz ? Maybe if we can extend before the code that reads and cuts, the script can take care of it.

[jrmerz]

Think we need to bring in @qjhart for this. The upper bound is either set to the current cost
until it's greater than the upper bound, then its set to the upper bound. The current code does not seem manipulate the cost curve, but I could be missing something.

[jdherman]

Notes from 9/23 meeting, we worked out the logic for how to handle these. Picture attached.
("slope" refers to any of the slopes for the full curve. LB/UB are the physical bounds for the whole link. Lowercase l/u are the lower/upper bounds on each piecewise piece in the matrix export).

We don't want to change the actual HOBBES data. This should all be happening in the matrix exporting script.

Edit: if you click the image, it will rotate properly..

Does this make sense?

[jrmerz]

@jdherman @msdogan Just a couple of clarification questions (for now).

1. Slope as defined by, k=0? Or if slope = 0 at k=0, check for k > 1
2. Top level 'else' statement includes 0 slope?
3. the expression 'for k = k'. Is that 'for all k'? or for last k? other?

[jdherman]

@jrmerz all good questions. Should have been clearer about these.

1. @msdogan mentioned that all "piecewise pieces" for a particular link should have slopes of the same sign. So checking the first one (k=0) should be enough ... but just for completeness, you may want to confirm that they're either ALL negative, or all positive.
2. I don't know how to handle 0 slope cases. For a negative slope link, this would only happen on the last piece -- and then you'd get an infinity if you tried to extend it to the x-intercept. So I'll hope that there just aren't any 0 slopes. Maybe flag it just to see? (Even 0 within tolerance would be a problem with the extension.
3. Sorry that should be k = K, 'the last k'.

[msdogan]

For the 2. question, we would want to do this extension only for links/nodes with piecewise cost. Most links/nodes don't have piecewise cost, so for them k=0, c=0 or constant c. I can think of two solutions:
1: Another if statement outside of "if cost is negative" saying if node/link has piecewise cost (something like k>0), then do the extension.
2: if elseif else. So, if slope is negative, then..., elseif slope is positive, then... else, break.

If we have piecewise representation (k>0) and c = 0, then we might have a problem.

[elliewhite-usgs]

@msdogan I don't think we will have the "piecewise representation (k>0) and c = 0" representation, because as Mimi said, a mixed linear-integer program will not be able to solve that. i.e. can't have the same slope on multiple pieces of the piecewise representation.

[jrmerz]

piecewise_processing_first_pass.txt

So I have implemented a first pass at this. I'm still trying to figure out how I will verify things are being done correctly. For now I am just adding console log to modifications made by the new peicewise logic.

I still need to review some of the code here, but take a glance a let me know how things look. This is for one Month, 1999-12-1.

[jrmerz]

Oops. Ignore that last file. Found issue. New file coming...

[jrmerz]

Always find the issue after you post... the first pass was not looking at the bounds correctly. Here is the updated output.

piecewise_processing_2nd_pass.txt

[msdogan]

@jrmerz this is great, but if possible can we get actual results (extension) from your script? We can look at examples for each case and plot them to see how they are looking.

[jrmerz]

Sure, let me know if you would like a different format. In the file, the costs array (in JSON format) is dumped after changes are made. Note, I realized the last file did not reflect your most resent calvin-network-data pull request. This current file does.

piecewise_processing_full_output.txt

[jrmerz]

@msdogan I have added the new version of the code. Couple notes, I have branched this code, so to test from the repository (ignore first step if you have already cloned repo).

git clone https://github.com/ucd-cws/calvin-network-tools
cd calvin-network-tools
git checkout extrapolated-cost
node bin/cnf.bin.js matrix --verbose --format=csv --start=1999-12 --stop=2000-02 --ts=. --to=network --max-ub=1000000 

Let me know if you have any issues.

[msdogan]

@jrmerz the code above is giving me an error

AD3+msdogan@CWS-Amargosa  /z/calvin-network-tools (extrapolated-cost)
$ node bin/cnf.bin.js matrix --verbose --format=csv --start=1999-12 --stop=2000-02 --ts=. --to=network --max-ub=1000000
*** CNF Arguments ***
  _allowUnknownOption: false
  _name: matrix
  _noHelp: false
  _description: Create a network matrix
  format: csv
  header: true
  start: 1999-12
  stop: 2000-02
  ts: .
  to: network
  maxUb: 1000000
  data: z:\calvin-network-data\data
  runtime: Z:\nodejs\node_modules\calvin-network-tools\HEC_Runtime
  verbose: true
[]
**********************

Running matrix command.

Z:\calvin-network-tools\nodejs\matrix\index.js:33
  hnf.split(config.data, {id:'prmname'}, config.nodes, function (subnet) {
      ^

TypeError: hnf.split is not a function
    at matrix (Z:\calvin-network-tools\nodejs\matrix\index.js:33:7)
    at module.exports.matrix (Z:\calvin-network-tools\nodejs\cmds\matrix.js:58:3
)
    at module.exports (Z:\calvin-network-tools\nodejs\run.js:41:45)
    at Z:\calvin-network-tools\nodejs\index.js:48:5
    at null._onTimeout (Z:\calvin-network-tools\nodejs\lib\checkVersion.js:15:5)

    at Timer.listOnTimeout (timers.js:92:15)

[jrmerz]

Oops. missed a step. npm install

So steps should be

git clone https://github.com/ucd-cws/calvin-network-tools
cd calvin-network-tools
git checkout extrapolated-cost
npm install
node bin/cnf.bin.js matrix --verbose --format=csv --start=1999-12 --stop=2000-02 --ts=. --to=network --max-ub=1000000 

[msdogan]

thanks @jrmerz it looks good in the first glance. I will take a look at cost curves for each case and let you know with my findings.

[msdogan]

Here are comparison of S09I05PD, S09I05 and pyvin.

Notes:

• S09I05PD does not have extended penalty data. It has original data sets
• S09I05 has HEC-PRM's penalty data extension (extrapolation to cover whole feasible range)
• pyvin has extended penalty data, but cost is represented as unit cost instead of paired data.
• S09I05PD and S09I05 were converted to matrix format (same as pyvin)

Findings:

• Results are looking good! Except for SR_SHA and SR_CLE, which both fall into k>0, slope<0, lb is not null category. For these two links, c is not printed in pyvin network export. Please see below for an example.

Examples for each "if" case:

This document shows examples for piecewise cost representation for each case discussed in https://github.com/ucd-cws/calvin-network-tools/issues/36#issuecomment-250806776

example format: origin_node terminal_node

Mustafa Dogan - 09/30/2016

if k > 0 # piecewise cost
    if slope < 0
        if LB != null
            example: SR_SHA.1999-12-31  SR_SHA.2000-01-31 , example: D5.1999-12-31  D73.1999-12-31
        else # LB = null
            example: C44.1999-12-31 C88.1999-12-31
        if UB != null # '!=' not equal to
            example: HU207.1999-12-31   CVPM08S.1999-12-31
        else # UB = null
            example: SR_NML.1999-12-31  D670.1999-12-31
    else # slope > 0
        if LB != null
            example: Should never happen!
        else # LB = null
            example: Should never happen!
        if UB != null # '!=' not equal to
            example: PMP_Banks.1999-12-31   D800.1999-12-31
        else # UB = null
            example: Should never happen!

example: D5.1999-12-31 D73.1999-12-31

k   S09I05PD            S09I05          pyvin       
    c   lb  ub  c   lb  ub  c   lb  ub
0   -4.11709    0   420.962 -4.117018244    0   420.960 -4.117085718    277 420.962
1   -4.11709    0   210.481 -4.117255796    0   210.480 -4.117086591    0   210.481
2   -4.11709    0   210.481 -4.116780692    0   210.480 -4.117085135    0   210.481
3   -4.11709    0   210.481 -4.117255796    0   210.480 -4.117085135    0   210.481

example: C44.1999-12-31 C88.1999-12-31

k   S09I05PD            S09I05          pyvin       
    c   lb  ub  c   lb  ub  c   lb  ub
0   -1.194332649    0   12.175  -1.194332649    0   12.175  -1.19436    0   12.175
1   -1.187612946    0   6.087   -1.187612946    0   6.087   -1.18768    0   6.087
2   -1.181205849    0   6.087   -1.181205849    0   6.087   -1.18115    0   6.087
3   -1.174141613    0   6.087   -1.174141613    0   6.087   -1.17393    0   4.191
4   -1.160315375    0   6.088   -1.160315375    0   6.088   -1.16050    0   0
5   -1.14654181 0   6.087   -1.14654181 0   6.087   -1.14651    0   0

example: HU207.1999-12-31 CVPM08S.1999-12-31

k   S09I05PD            S09I05          pyvin       
    c   lb  ub  c   lb  ub  c   lb  ub
0   -1265.000   0   0.010   -1265.000   0   0.080   -1265.000   0   0.010
1   -1263.000   0   0.010   -1263.000   0   0.010   -1263.000   0   0.010
2   -663.000    0   0.010   -663.000    0   0.010   -663.000    0   0.010
3   -235.500    0   0.010   -235.000    0   0.010   -235.500    0   0.010
4   -233.500    0   0.010   -234.000    0   0.010   -233.500    0   0.010
5   -52.000 0   0.010   -52.000 0   0.010   -52.000 0   0.080
6               0   0   999999999.870           

SR_NML.1999-12-31 D670.1999-12-31

k   S09I05PD            S09I05          pyvin       
    c   lb  ub  c   lb  ub  c   lb  ub
0   -0.91320    0   177.080 -0.91320    0   177.080 -0.913166155    0   177.084
1   -0.90796    0   88.550  -0.90796    0   88.550  -0.908059384    0   88.542
2   -0.90309    0   88.540  -0.90309    0   88.540  -0.903072947    0   88.542
3   -0.89756    0   88.540  -0.89756    0   88.540  -0.897552903    0   88.542
4   -0.88740    0   88.540  -0.88740    0   88.540  -0.887282851    0   88.542
5   -0.87645    0   88.550  -0.87645    0   88.550  -0.876587022    0   88.542
6               0   0   999999380.200           

PMP_Banks.1999-12-31 D800.1999-12-31

k   S09I05PD            S09I05          pyvin       
    c   lb  ub  c   lb  ub  c   lb  ub
0   16.46473    0   47.210  16.46418435 0   47.214  16.46557933 0   47.210
1   16.62219    0   47.220  16.62345914 0   47.214  16.62198221 0   47.220
2   16.78034    0   47.210  16.7796323  0   47.212  16.78034315 0   47.210
3   16.87421    0   47.220  16.87420584 0   47.220  16.87335875 0   47.220
4   16.94556    0   47.210  16.94556238 0   47.210  16.9457742  0   47.210
5   17.01334    0   47.210  17.01334463 0   47.210  17.01419191 0   47.210
6   17.06480    0   47.220  17.06480305 0   47.220  17.06480305 0   47.220
7   17.10655    0   47.210  17.10654522 0   47.210  17.10569795 0   47.210
8   17.14104    0   47.220  17.14104193 0   47.220  17.14125371 0   47.220
9   17.21034    0   47.210  17.20807438 0   999999575.070   17.20948951 0   5.886

PROBLEM
SR_SHA.1999-12-31 SR_SHA.2000-01-31

SR_SHA.1999-12-31,SR_SHA.2000-01-31,0,,0.999061850013677,509.89398200000005,509.89398200000005
SR_SHA.1999-12-31,SR_SHA.2000-01-31,1,,0.999061850013677,729.0797729999999,729.0797729999999

[msdogan]

@jrmerz could you please let me know if there is any link falling into "should never happen!" category? Thanks

[jrmerz]

@msdogan The following links are in the 'should never happen' category:

D550-PMP_CC1: s>0 && !UB, Extending k to ub=1000000000
 WARNING This should never happen!
D550-PMP_CC1: s>0 && !UB, Extending k to ub=1000000000
 WARNING This should never happen!
C309-PMP_Old_R: s>0 && !UB, Extending k to ub=1000000000
 WARNING This should never happen!
C309-PMP_Old_R: s>0 && !UB, Extending k to ub=1000000000
 WARNING This should never happen!
D528-PMP_Mallar: s>0 && !UB, Extending k to ub=1000000000
 WARNING This should never happen!
D528-PMP_Mallar: s>0 && !UB, Extending k to ub=1000000000
 WARNING This should never happen!

[msdogan]

okay, well. Sometimes it can happen :) Those are outliars. Here what metada says for penalty on C309-PMP_Old_R

"CCWD has three sources of water, all of which have monthly varying water quality. The penalty set is to reflect the salinity damage of each source (based on TDS above the 250 mg/l at 50 cents/ mg/l TDS) plus the fixed treatement and local distrubtion and conveyance costs ($299/af). REVISED 12-27-2002 MJ - removed $299 from the paired data cost, so it reflects only salinity damage; aslo recomputed see revised file on Roselyn in folder "G03G07""

I compared your extension to HEC-PRM, and it complies with that. So, we are good. You can remove the warning for this case :)

Only thing is that UB should be max_UB without adding a new link, because there is no physical upper bound. So, in the example below, there are two links with the same cost. In this case, the solver wouldn't know which link to use.

C309.1999-12-31 PMP_Old_R.1999-12-31    0   81.49999811 1   0   15.38
C309.1999-12-31 PMP_Old_R.1999-12-31    1   81.49999811 1   0   1000000
C309.2000-01-31 PMP_Old_R.2000-01-31    0   51.00000033 1   0   15.38
C309.2000-01-31 PMP_Old_R.2000-01-31    1   51.00000033 1   0   1000000

[jrmerz]

k, pushed to npm @ calvin-network-tools@1.0.48

[msdogan]

Closing the issue. The solution from pyvin is infeasible but we are going to tackle with that in a separate issue.

[jrmerz]

Ok, latest code pushed to npm calvin-network-tools@1.0.49

[msdogan]

@jrmerz @jdherman reopening this issue. I think it is related to piecewise cost representation for the first piece. I run a one-year model between Oct-2002 and Sep-2003 in debug mode and it was feasible. I created another dataset for Oct-1998 and Sep-2003 but pyomo gave an error because for some months lb > ub. Below is the list, it is happening for SR_SHA, SR_CLE

SR_SHA.2001-01-31   SR_SHA.2001-02-28   0   -6.82366224 0.999185861 2000    532.688904
SR_SHA.2001-02-28   SR_SHA.2001-03-31   0   -5.712019564    0.998986365 2000    634.428528
SR_SHA.2001-03-31   SR_SHA.2001-04-30   0   -5.452069302    0.997972736 2000    736.177979
SR_SHA.2001-04-30   SR_SHA.2001-05-31   0   -4.848080055    0.99754679  2000    759.528381
SR_SHA.2001-05-31   SR_SHA.2001-06-30   0   -4.874155583    0.994807825 2000    762.174988
SR_SHA.2001-06-30   SR_SHA.2001-07-31   0   -4.85020368 0.994780871 2000    762.172607
SR_SHA.2001-07-31   SR_SHA.2001-08-31   0   -6.454045852    0.993923596 2000    732.284729
SR_SHA.2001-08-31   SR_SHA.2001-09-30   0   -6.339085112    0.994257875 2000    733.677673
SR_SHA.2001-09-30   SR_SHA.2001-10-31   0   -6.652135027    0.9957298   2000    561.510498
SR_SHA.2001-10-31   SR_SHA.2001-11-30   0   -7.003443177    0.996738046 2000    514.402527
SR_SHA.2001-11-30   SR_SHA.2001-12-31   0   -6.971974325    0.998474161 2000    493.441956
SR_SHA.2001-12-31   SR_SHA.2002-01-31   0   -7.316334977    0.99906185  2000    509.893982
SR_SHA.2002-01-31   SR_SHA.2002-02-28   0   -6.82366224 0.99919125  2000    532.688904
SR_SHA.2002-02-28   SR_SHA.2002-03-31   0   -5.712019564    0.99887314  2000    634.428528
SR_SHA.2002-03-31   SR_SHA.2002-04-30   0   -5.452069302    0.998333975 2000    736.177979
SR_SHA.2002-04-30   SR_SHA.2002-05-31   0   -4.848080055    0.997131636 2000    759.528381
SR_SHA.2002-05-31   SR_SHA.2002-06-30   0   -4.874155583    0.996101825 2000    762.174988
SR_SHA.2002-06-30   SR_SHA.2002-07-31   0   -4.85020368 0.99485635  2000    762.172607
SR_SHA.2002-07-31   SR_SHA.2002-08-31   0   -6.454045852    0.993670186 2000    732.284729
SR_SHA.2002-08-31   SR_SHA.2002-09-30   0   -6.339085112    0.994462761 2000    733.677673
SR_SHA.2002-09-30   SR_SHA.2002-10-31   0   -6.652135027    0.995562661 2000    561.510498
SR_SHA.2002-10-31   SR_SHA.2002-11-30   0   -7.003443177    0.99744435  2000    514.402527
SR_SHA.2002-11-30   SR_SHA.2002-12-31   0   -6.971974325    0.998587386 2000    493.441956
SR_SHA.2002-12-31   SR_SHA.2003-01-31   0   -7.316334977    0.999024111 2000    509.893982
SR_CLE.2001-01-31   SR_CLE.2001-02-28   0   -1.642629326    0.999780165 1100    392.994141
SR_CLE.2001-02-28   SR_CLE.2001-03-31   0   -1.462461308    0.999590998 1100    392.844177
SR_CLE.2001-03-31   SR_CLE.2001-04-30   0   -1.520100798    0.998456023 1100    375.905457
SR_CLE.2001-04-30   SR_CLE.2001-05-31   0   -1.376631074    0.997960113 1100    377.476318
SR_CLE.2001-05-31   SR_CLE.2001-06-30   0   -1.388836865    0.995869102 1100    374.334656
SR_CLE.2001-06-30   SR_CLE.2001-07-31   0   -1.369219249    0.995562352 1100    392.877136
SR_CLE.2001-07-31   SR_CLE.2001-08-31   0   -1.673263732    0.994892613 1100    522.918518
SR_CLE.2001-08-31   SR_CLE.2001-09-30   0   -1.631388432    0.995199363 1100    538.842712
SR_CLE.2001-09-30   SR_CLE.2001-10-31   0   -1.523674962    0.99674845  1100    541.927307
SR_CLE.2001-10-31   SR_CLE.2001-11-30   0   -1.564700716    0.998185065 1100    537.30011
SR_CLE.2001-11-30   SR_CLE.2001-12-31   0   -1.671675597    0.999473415 1100    382.097778
SR_CLE.2001-12-31   SR_CLE.2002-01-31   0   -1.753488658    0.999764827 1100    392.415466
SR_CLE.2002-01-31   SR_CLE.2002-02-28   0   -1.642629326    0.999790385 1100    392.994141
SR_CLE.2002-02-28   SR_CLE.2002-03-31   0   -1.462461308    0.9995092   1100    392.844177
SR_CLE.2002-03-31   SR_CLE.2002-04-30   0   -1.520100798    0.998936602 1100    375.905457
SR_CLE.2002-04-30   SR_CLE.2002-05-31   0   -1.376631074    0.997627798 1100    377.476318
SR_CLE.2002-05-31   SR_CLE.2002-06-30   0   -1.388836865    0.996508167 1100    374.334656
SR_CLE.2002-06-30   SR_CLE.2002-07-31   0   -1.369219249    0.995593027 1100    392.877136
SR_CLE.2002-07-31   SR_CLE.2002-08-31   0   -1.673263732    0.99478525  1100    522.918518
SR_CLE.2002-08-31   SR_CLE.2002-09-30   0   -1.631388432    0.995347627 1100    538.842712
SR_CLE.2002-09-30   SR_CLE.2002-10-31   0   -1.523674962    0.99668199  1100    541.927307
SR_CLE.2002-10-31   SR_CLE.2002-11-30   0   -1.564700716    0.998420238 1100    537.30011
SR_CLE.2002-11-30   SR_CLE.2002-12-31   0   -1.671675597    0.99950409  1100    382.097778
SR_CLE.2002-12-31   SR_CLE.2003-01-31   0   -1.753488658    0.99974949  1100    392.415466

[msdogan]

the list above was only when lb > ub.
Here is some months from 1-year run for SR_SHA. For example for SR_SHA.2002-10-31 SR_SHA.2002-11-30 above output is wrong but below (1-year) is correct. Not sure how it is happening.

SR_SHA.2002-10-31   SR_SHA.2002-11-30   0   -7.003443177097767  0.9974443500372581  514.402527  514.402527
SR_SHA.2002-10-31   SR_SHA.2002-11-30   1   -2.973775927491251  0.9974443500372581  737.3889770000001   737.3889770000001
SR_SHA.2002-10-31   SR_SHA.2002-11-30   2   -1.4663235690271152 0.9974443500372581  448.20849599999997  2148.2084959999997
SR_SHA.2002-11-30   SR_SHA.2002-12-31   0   -6.971974324777523  0.9985873855105941  493.441956  493.441956
SR_SHA.2002-11-30   SR_SHA.2002-12-31   1   -3.0555296683697994 0.9985873855105941  700.8379510000001   700.8379510000001
SR_SHA.2002-11-30   SR_SHA.2002-12-31   2   -1.4789632953596268 0.9985873855105941  505.7200929999998   2057.720093
SR_SHA.2002-12-31   SR_SHA.2003-01-31   0   -7.316334976865838  0.9990241105042272  509.89398200000005  509.89398200000005
SR_SHA.2002-12-31   SR_SHA.2003-01-31   1   -3.127525157662001  0.9990241105042272  729.0797729999999   729.0797729999999
SR_SHA.2002-12-31   SR_SHA.2003-01-31   2   -1.5365093619413377 0.9990241105042272  461.026245  2129.026245
SR_SHA.2003-01-31   SR_SHA.2003-02-28   0   -6.823662240203149  0.9991319395226552  532.688904  532.688904
SR_SHA.2003-01-31   SR_SHA.2003-02-28   1   -2.666343546632485  0.9991319395226552  889.55487   889.55487
SR_SHA.2003-01-31   SR_SHA.2003-02-28   2   -1.2979171742850106 0.9991319395226552  577.7562259999999   2405.756226
SR_SHA.2003-02-28   SR_SHA.2003-03-31   0   -5.712019564164365  0.9989432355054063  634.428528  634.428528
SR_SHA.2003-02-28   SR_SHA.2003-03-31   1   -1.8320108363125103 0.9989432355054063  1324.7074579999999  1324.7074579999999
SR_SHA.2003-02-28   SR_SHA.2003-03-31   2   -1.0003759700998833 0.9989432355054063  40.864014000000225  2082.864014

[msdogan]

the lb - ub issue I just posted above a few days ago also related to piecewise cost representation but not extension. Let's solve extension issue first, and then we can try to solve lb - ub issue for SR_SHA and SR_CLE.
Once again, below figure shows how HEC-PRM extrapolates the cost functions (piecewise linear penalties).
So, the issue here is when the slope is negative and there is not upper bound, we try to intersect x-axis by extending the last piece's ub. This is the case where slope is negative, and ub does not exist. However, HEC-PRM does not try to intersect x-axis. If you look at the figure below, in (a), (b), (c) and (d), HEC-PRM adds new pieces with zero slope. It does not try to cross x-axis.

In our CALVIN/pyvin system this becomes an issue for reservoirs releases where there is hydropower. In piecewise cost, there is maximum capacity that represents turbine capacity, at which cost is zero. However, we do not put a upper bound on reservoir releases. The reason may be that turbine is not only release, there is also spilling, which I don't think possible to represent in linear programming, at least we don't represent right now.

When slope is positive, there is no issue, our extension complies with (e) and (f).
Suggestions for right hand side extension when slope is negative (a, b, c, d):

• If physical ub exists, add a new piece with a cost of zero and ub for this new piece is physical ub minus ub on last the last piece. If not, add a new piece with ub is max ub.
• The second suggestion is the same as first one, but first extend to intersect x-axis first, and then extend to physical ub or max ub

First suggestion is what HEC-PRM is doing even if we add a piece with zero cost. Since other pieces will have negative cost, it will use those first before using arc with zero cost. We have to make sure to put physical upper bound to prevent flooding or dumping water (zero cost) for necessary links, which is why we are adding target deliveries as upper bounds for ag and urban.
The second suggestion requires a bit more thinking because there will be more sub-cases with more if statements

@jdherman @mimijenkins1 @josue-medellin Please let us know what you are thinking and then we can proceed. Thanks.

[jdherman]

Thanks for summarizing @msdogan .

If I'm understanding right, it sounds like this problem is specific to reservoirs (SR_*), because reservoirs have an upper bound in the dataset that does NOT reflect the physical release capacity, only the turbine capacity. Our logic before did not account for this.

So when we enforce this upper bound, the reservoir can't release enough water to meet mass balance constraints.

If this is the case, then the only modification in the matrix export script is to not enforce ub for reservoir release links.

My guess is this also is causing the lb > ub issue for SHA and CLE -- this would occur for example if some required environmental releases (lb) exceed the turbine capacity (misinterpreted as ub).

@jrmerz does this make sense? This would be:

if slope negative
  if UB exists
    if link not reservoir
      enforce u = UB for last k
...

[msdogan]

@jdherman this is a issue for links where there is a downward sloping cost (negative) but there is no physical upper bound. Reservoir release is one example. And it makes sense not to put a physical upper bound because you can release through turbines or spill with a zero benefit. And usually spilling capacities are big. As you said, when we enforce turbine capacity for reservoir releases (from piecewise curve), it does not release enough water.
Assuming this is only happening for reservoirs if the link is not reservoir then enforce u=UB. What if link is a reservoir? :) We need to something for reservoir release links. Right now, we know this is for reservoirs but it might be happening in other parts of the network. I think making it general would be better.

[jdherman]

Ok -- I'm fine with making it general. But we still have the problem of identifying where this is happening, because we have to somehow separate out the links for which the upper bound value in the dataset is not representative of the real-world upper bound. How to do that? It must be based on a specific type of nodes/links.

[jrmerz]

Could this be a declarative flag on the data? So Any link/node can declaratively specify 'representative' v 'real-world' upper bound? And then code would look like:

if slope negative
  if UB exists and bound not representative
      enforce u = UB for last k

This would assume real world by default. So the flag would be representative = true.

[jdherman]

Something like this, yea. Our problem in this case is that the "upper bound" in the data for reservoirs actually represents the turbine capacity, not the flow capacity constraint. So it is representative of something, just not what we want it to be.

@msdogan can you think of other types of links where this might be happening? Or is it only reservoir release links?

[msdogan]

I can't think of any other types of links but @jrmerz could you do a formal check for this case: piecewise cost, slope negative and no physical upper bound?

[jrmerz]

Here is the current list meeting the above conditions:

PWP_Mojave-SR_SLW, SR_NML-D670, SR_LL_ENR-SR_DNP, SR_HTH-C44, SR_DNP-D662, SR_MCR-D642, D541-Surp_Delta, SR_SHA-D5, SR_ORO-C23, SR_ENG-C28, SR_BUL-C27, D5-D73, SR_TAB-C25, SR_FOL-D9, D9-D85, C23-C25, SR_CLE-D94, PWP_AAC-C151

[jdherman]

Current thought:

if slope negative:
  if UB exists:
    add new piece up to UB at cost 0 (flat slope)
  else:
    add new piece up to max_UB at cost 0 (flat slope)

[jdherman]

Copying this from Mustafa's email just so we have a record of it:

Okay, while checking the network, I found another bug in big "if statement" logic. https://cloud.githubusercontent.com/assets/1936620/18797140/505be90c-8181-11e6-84f6-d0ef0ae7e0eb.JPG

For left-hand side extension for slope negative and physical lb does not exist case, our statement says: enforce lb=0 for k=0 case. But does not say anything for the ub. Before the extension, ub on the k=0 case (sorry I forgot start from 0 index in the sketch below) is b - a. After extension it should be b, but it is still b - a in network matrix because we forgot to add this in the if statement.

Sorry Justin, but could you update your code again to include that. Hopefully this would be the last modification for cost extension.

[jrmerz]

Logic updated. But the cost of slope negative, Physical Lower Bound does not exist, we only modify lb if lb wasn't already equal to 0. Currently, no nodes/links meet this condition. In fact no nodes/links have slope negative and an undefined Physical Lower Bound.

[lzhangisu2011]

Just saw another example might be related to piecewise cost representation issue:

Infeasibility is for reservior node SR_MCR in the final month.

Here is all inflows to this node:
i j k cost amplitude lower_bound upper_bound
SR_MCR.2002-11-30 SR_MCR.2002-12-31 0 -2.309251768 0.999606565 115 126.214554
SR_MCR.2002-11-30 SR_MCR.2002-12-31 1 -1.473401553 0.999606565 0 182.383865
SR_MCR.2002-11-30 SR_MCR.2002-12-31 2 -1.009400653 0.999606565 0 254.401581
INFLOW.2002-12-31 SR_MCR.2002-12-31 0 0 1 33.66198 33.66198

The outflow of this node is:
i j k cost amplitude lower_bound upper_bound
SR_MCR.2002-11-30 SR_MCR.2002-12-31 2 -1.009400653 0.999606565 0 254.401581
SR_MCR.2002-12-31 FINAL 0 0 1 633.792908 633.792908

There is about 37.1309 difference between inflow and outflow.

[msdogan]

I found the bug, the upper bound on the first piece (k=0) of SR_MCR storage link should be 242.287 not 126.214554.

SR_MCR.2002-11-30 SR_MCR.2002-12-31 0 -2.309251768 0.999606565 115 126.214554

We probably missed this in the big if statement. I am going to the class now but I will try to find where to fix once I get a chance. Thanks @lzhangisu2011 for finding it.

[msdogan]

@jrmerz I am closing this issue. It is getting too large and hard to follow. We will open a new one if we get any problem with piecewise cost representation in the future.

[jrmerz]

agreed.