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

[msdogan]

@jdherman @jrmerz just to clarify one thing. When we say UB exists in the if statement, we mean physical UB, right? Not UB on the link on penalty curve.

[jrmerz]

Have to show my ignorance again, at least this repo isn't public :)

So LB & UB means the bound property on a node/link, as defined by a file such as: https://github.com/ucd-cws/calvin-network-data/blob/master/data/c88-c78/UBM.csv

[jrmerz]

... at least that's whats being checked in the code

[jdherman]

Yes, that's what I was thinking (physical UB on the link). But if there is no physical UB for these links, then where are those UB's coming from (the 15.38 value)? Target deliveries? Should that still count as a UB? I don't know anymore.

I agree that the extension should happen to existing links without creating new links.

[msdogan]

@jdherman 15.38 is from penalty curve, but there is no physical upper bound on this link.

@jrmerz yes that example is physical upper bound on C88-C78.

[jrmerz]

@msdogan I have pushed fixes for the positive slope, no upper bound issue. I have also pushed a fix for nodes that lack flow but still have inflow. Pushed to 'extrapolated-cost' branch. Please pull the 'extrapolated-cost' branch and test. Once verified I will merge extrapolated-cost back to master and publish changes to NPM.

[msdogan]

@jrmerz fixes for D550-PMP_CC1, C309-PMP_Old_R and D528-PMP_Mallar are looking good.

For #37, those inflows are still missing in matrix export.

One last thing, this is probably not included in this push, but cost (c) values are still missing for SR_SHA and SR_CLE

SR_SHA.1999-12-31   SR_SHA.2000-01-31   0       0.99906185  509.893982  509.893982
SR_SHA.1999-12-31   SR_SHA.2000-01-31   1       0.99906185  729.079773  729.079773

SR_CLE.1999-12-31   SR_CLE.2000-01-31   0       0.999764827 392.415466  392.415466
SR_CLE.1999-12-31   SR_CLE.2000-01-31   1       0.999764827 308.909302  308.909302

[jrmerz]

@msdogan k=2 for SR_SHA is being left out because the current upper bound (cub in code) is 0. And this is explicitly ignored for storage. https://github.com/ucd-cws/calvin-network-tools/blob/master/nodejs/matrix/storage.js#L103

It looks like issue #32 attacked this same condition for links by removing this condition: https://github.com/ucd-cws/calvin-network-tools/blob/master/nodejs/matrix/link.js#L98

Should I do the same for storage?

[msdogan]

@jrmerz there is a few issues with these two storage nodes: 1- cost values are not written in the matrix export 2- lower bound and upper bound are the same, which shouldn't be. 3- Finally, I wonder why the code says cub = 0.

For SR_SHA and SR_CLE, there are lower bound, upper bound constraints and storage cost. In our big `if statement', this is the case where slope is negative, lb exists, and ub exists.

[jrmerz]

@msdogan ok, found the part of the issue with the cost. It was accessing the wrong variable, this should have actually been costs[k].cost, will push change in a second. As for the calculations for cub. This is not code I am super familiar with. I have cleaned up the code a little bit here: https://github.com/ucd-cws/calvin-network-tools/blob/master/nodejs/matrix/storage.js#L90

And, here is an attempt to show how SR_SHA storage is being calculated (note clb and cub are the variables added to the matrix):

Starting stepBounds: {"LB":2000,"UB":3368}
k=0
  costs[0] = {"cost":-7.316334976865838,"lb":2000,"ub":509.89398200000005}
  (costs[0].ub || 0) <= stepBounds.LB
    Setting clb = costs[0].ub or 0
  Reducing stepBounds.LB by clb.
  costs[0].ub !== null && costs[0].ub <= stepBounds.LB
    Setting cub = costs[0].ub
  Reducing stepBounds.UB by cub.

  Step Values:
    k = 0
    cost = -7.316334976865838
    clb = 509.89398200000005
    cub = 509.89398200000005
    stepBounds.UB = 2858.106018
    stepBounds.LB = 1490.106018
  row added

k=1
  costs[1] = {"cost":-3.127525157662001,"lb":0,"ub":729.0797729999999}
  (costs[1].ub || 0) <= stepBounds.LB
    Setting clb = costs[1].ub or 0
  Reducing stepBounds.LB by clb.
  costs[1].ub !== null && costs[1].ub <= stepBounds.LB
    Setting cub = costs[1].ub
  Reducing stepBounds.UB by cub.

  Step Values:
    k = 1
    cost = -3.127525157662001
    clb = 729.0797729999999
    cub = 729.0797729999999
    stepBounds.UB = 2129.026245
    stepBounds.LB = 761.026245
  row added

k=2
  costs[2] = {"cost":-1.5365093619413377,"lb":0,"ub":3368}
  (costs[2].ub || 0) > stepBounds.LB && costs[2].lb <= stepBounds.LB
    Setting clb = stepBounds.LB
  Reducing stepBounds.LB by clb.
  costs[2].ub === null || costs[2].ub > stepBounds.LB
    Setting cub = stepBounds.LB
  Reducing stepBounds.UB by cub.

  Step Values:
    k = 2
    cost = -1.5365093619413377
    clb = 761.026245
    cub = 0
    stepBounds.UB = 2129.026245
    stepBounds.LB = 0
  row not added, cub <= 0

[jrmerz]

Also note, the above code happens after the logic described earlier in this issues.

[msdogan]

@jrmerz okay, the code seems to be doing good. I understand why it constraints. Physical lower bound is 2000, so it constraints all pieces until 2000. But, it does not show anything between 2000 and 3368, which is physical upper bound.

k=0 and k=1 seems correct. I think problem here is that 2000 lb corresponds 3rd piece, k=2. We may need to modify the code. For k=2, we want to constraint from end of k=1 to 2000 lb, but after that there will be unconstained (I mean constrained but ub > lb, there will a range with a cost).

In the figure below vertical continuous lines show lb and ub, and k=0,1,2 pieces are between dashed lines.

@jdherman @mimijenkins1 I don't think we thought about this. So, what do you think? :) Do you think we need to add a new piece when a penalty line is crossed with lower bound? Or do we enforce lower bound on third piece and don't put first and second pieces in network matrix?

[jdherman]

Yea, we missed this case in our logic.

I can see two equally good options for doing this (and neither one should require an extra piece):

1. Set the first two pieces as fixed flows (lb = ub) because you know they must be full. Then the lb can be specified on the third piece. Or,
2. Don't include the first two pieces in the matrix at all, and set the lower bound on the third piece as the physical lower bound.

Either one is fine by me, but it might be clearer to do (1) and keep all the pieces in the matrix just so we can keep track of them. Since they're all full (X=ub) there should be no effect on shortage cost.

[jrmerz]

UCSC group call on 10/19 has decided to try 1) above.

[msdogan]

@jrmerz in the last fix (https://github.com/ucd-cws/calvin-network-tools/commit/24a04d037c98c544467b036c9e710f8d0bb6d0ac), last piece is not in matrix export

SR_SHA.1999-12-31   SR_SHA.2000-01-31   0       0.99906185  509.893982  509.893982
SR_SHA.1999-12-31   SR_SHA.2000-01-31   1       0.99906185  729.079773  729.079773

based on your code above:

k=2
  costs[2] = {"cost":-1.5365093619413377,"lb":0,"ub":3368}
  (costs[2].ub || 0) > stepBounds.LB && costs[2].lb <= stepBounds.LB
    Setting clb = stepBounds.LB
  Reducing stepBounds.LB by clb.
  costs[2].ub === null || costs[2].ub > stepBounds.LB
    Setting cub = costs[2].ub - costs[1].ub - costs[0].ub
  Reducing stepBounds.UB by cub.

[msdogan]

This was the only way I can think of. There might be some other clever ways.

[msdogan]

Thanks @jrmerz it seems to work now

[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.