Include Transport Tutorial

[glatterf42]

The transport tutorial seems to be a standard tutorial for showcasing usage of solvers and related packages. GAMS has it as well as pyomo and, maybe in particular, ixmp has its own version of it. We are currently rewriting ixmp to be all-python and we are thinking of moving away from GAMS in the near future. Therefore, we are considering how our implementation choices now might enable or harm a future move to solver interfaces such as linopy.

After reading the documentation and missing the transport tutorial for comparison, I decided to try adding it myself, which is what this PR does. While I followed your contributing guidelines for the most part, I am not sure where in the docs this tutorial should be linked, so I did not complete step 4 yet. Please also check that:

• this is actually the best way to model the transport problem in linopy (coordinated with @aurelije)
• the text I have written is correct (includes feedback)
• keeping the GAMS/pyomo commands for comparison is not too confusing (restructured based on feedback for readability)

Regarding this last point, I should explain that most comments follow this structure: lines on top are usually from GAMS and marked with single #, lines below that (marked with ##) are from pyomo. If additional lines appear using just one # again, they are new comments intended to aid with the tutorial.

I have one note regarding your contributing guidelines: the pre-commit hook you linked did not work out of the box for me (using Ubuntu 22.04) because bash is not the default shell there and needs to be specified.

Other than that, I want to commend your documentation. In general, it was not a big problem to write up this tutorial without prior experience with pyomo or linopy. The most confusing part for me was the lack of sets and parameters as model attributes, it might be worth expanding your docs (migrating from pyomo) in that regard. Or do you plan to include sets and parameters as model attributes in the future?

[codecov]

Codecov Report

Patch coverage has no change and project coverage change: -0.05% ⚠️

Comparison is base (d9abfec) 87.03% compared to head (4612b7b) 86.98%.
Report is 30 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #149      +/-   ##
==========================================
- Coverage   87.03%   86.98%   -0.05%     
==========================================
  Files          14       14              
  Lines        3101     3120      +19     
  Branches      707      711       +4     
==========================================
+ Hits         2699     2714      +15     
- Misses        294      296       +2     
- Partials      108      110       +2     

see 4 files with indirect coverage changes

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[FabianHofmann]

@glatterf42 thank you very much for the PR! I will have a look at it asap, but I cannot say when that is (likely end of the week)

[aurelije]

@glatterf42 thanks for this addition. I myself wanted to to exactly the same model which is "Hello World" of OR. I have implemented it on my pc but didn't created nice explanations as you have.

I have took a look so here are my impressions. I don't like usage of numpy. Yes it is a working horse in ML and well known to all people in DS, it is not so hard to learn it and understand the code even if you have never worked with it. But it is just one more cognitive load, one more library to include and it brings nothing. Someone who is only OR/Math optimization oriented (has no ML/DS background) and maybe struggles with Python itself (is not developer/programmer) will scratch it's head.

As a Software engineer I also do not like usage of internal things like NumPy (on top of which xarray is built). There is a principle that you should talk only to your immediate neighbor, but not to neighbor of neighbor. All things done with numpy could be done using xarray at the first place.

For some lines I see over-complication. Maybe that is just a matter of stile (explicit vs implicit) but for example sum called on objective function is not needed at all, linopy would enforce sum by itself.

Here is what I had as a implementation, I am not saying it is better but maybe I can get some review to learn linopy better and maybe you can get some inspiration for doing this example a bit different. I was more respective to clean code implementation and python/linopy possibilities showcase than to mimic GAMS example I was reimplementing. So for comments, names of the datasets, names of the columns... I tried to preserve connection to GAMS but for names of python variables I have used more descriptive naming than COBOL/Fortran/Maths inspired one letter names like i, j, f...):

from linopy import Model
import xarray as xr

"""This problem finds a least cost shipping schedule that meets
requirements at markets and supplies at factories.


Dantzig, G B, Chapter 3.3. In Linear Programming and Extensions.
Princeton University Press, Princeton, New Jersey, 1963.

This formulation is described in detail in:
Rosenthal, R E, Chapter 2: A GAMS Tutorial. In GAMS: A User's Guide.
The Scientific Press, Redwood City, California, 1988.

The line numbers will not match those in the book because of these
comments.

Keywords: linear programming, transportation problem, scheduling"""

UNIT_PRICE_PER_K_MILES = 90
plants = {"canning plants": ['seattle', 'san-diego']}
markets = {"markets": ['new-york', 'chicago', 'topeka']}
capacity = xr.DataArray([350, 600], coords=plants, name="capacity of plant i in cases")
demand = xr.DataArray([325, 300, 275], coords=markets, name="demand at market j in cases")
unit_cost_in_k_dollars_per_transport_group = xr.DataArray(
    [[2.5, 1.7, 1.8], [2.5, 1.8, 1.4]],
    coords=plants | markets,
    name='transport cost in thousands of dollars per case'
) * (UNIT_PRICE_PER_K_MILES / 1000)


m = Model(force_dim_names=True)

transport = m.add_variables(lower=0,
                            coords=unit_cost_in_k_dollars_per_transport_group.coords,
                            name="shipment quantities in cases")

print("\n\ntransport vars:")
print(transport)

m.add_objective(unit_cost_in_k_dollars_per_transport_group * transport)

respect_capacity_of_production = transport.sum(dims='markets') <= capacity
m.add_constraints(respect_capacity_of_production, name='observe supply limit at plant i')

respect_market_demand = transport.sum(dims='canning plants') >= demand
m.add_constraints(respect_market_demand, name='satisfy demand at market j')
print("\n\nmodel before solving:")
print(m)

#m.solve(solver_name='gurobi', io_api='direct')
m.solve(solver_name='cbc')
print("\n\nmodel solved:")
print(m)
print("\n\nmodel solution:")
print(m.solution)
sol = m.solution.to_dataframe()
print("\n\nsolution dataframe:")
print(sol)

[FabianHofmann]

Thanks @aurelije, for the comment. I think the code you provided is very clean and nice. So I would definitely make sense to include it.

[glatterf42]

Thanks for the comment, @aurelije, and sorry for taking so long, everyone. I've now updated the tutorial according to your suggestions, e.g. getting rid of numpy etc. However, I kept the variable names analogous to GAMS/Pyomo to keep comparison with these tutorials easy.
There is one question that maybe @FabianHofmann could clarify: where in the docs can I find that .sum() is default for the objective? It might be nice to add that to explain why defining obj = c * x is enough in this case.

[FabianHofmann]

Thanks for the comment, @aurelije, and sorry for taking so long, everyone. I've now updated the tutorial according to your suggestions, e.g. getting rid of numpy etc. However, I kept the variable names analogous to GAMS/Pyomo to keep comparison with these tutorials easy. There is one question that maybe @FabianHofmann could clarify: where in the docs can I find that .sum() is default for the objective? It might be nice to add that to explain why defining obj = c * x is enough in this case.

Hey @glatterf42, that is correct. Since linopy assumes only one objective function, it automatically takes the full sum. So, either making a small note or explicitly taking the sum should be fine :)

I really like the example, great work! Would it be possible to pull some of the comments in the code cells into the text? For the comparison with the pyomo code you could insert boxes like

Equivalent in Pyomo

model.i = Set(initialize=['seattle','san-diego'], doc='Canning plans')
model.j = Set(initialize=['new-york','chicago', 'topeka'], doc='Markets')

[glatterf42]

Thanks for the clarification and suggestion, let me know what you think of the new structure :)

[FabianHofmann]

Nice work! Just roughly scanning, I would be happy to merge it. Let me know when you think it is ready.

[glatterf42]

Thanks! From my side, you are free to merge whenever you choose :)

[FabianHofmann]

Thanks @glatterf42 for this nice contribution! I will make sure it appears on the doc asap.