Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

bugfixes for more robust feasibility #744

Merged
merged 32 commits into from
Nov 11, 2024
Merged

Conversation

flohump
Copy link
Contributor

@flohump flohump commented Nov 5, 2024

🐦 Description of this PR 🐦

This PR aims to fix infeasibilies that I have repeatedly encountered for single runs when testing the current develop with the full set of FSEC runs and runs with H16 resolution.
The fix in 10_land for the land transtion matrix is to use only variables (vm_land) and variables levels (vm_land.l), and not a mix of variables and parameters. The issue is that GAMS seems to have different accuracies for variables and parameters. However, vm_land.l is updated between time steps based on parameters. Therefore, the balance variables are still needed.

This PR has been successfully tested with test_runs.R, the full set of FSEC runs and some H16EU runs with rev4.114. Tests with rev4.115 are underway.

  • 10_land bugfix land transition matrix for improved feasibility (variables and parameters have different accuracy)
  • 44_biodiversity avoid division by zero and improved consistency between realisations for fixing variables
  • 70_livestock bugfix scaling.gms file in wrong folder
  • modules update of scaling factors in several modules

🔧 Checklist for PR creator 🔧

  • Label pull request from the label list.

    • Low risk: Simple bugfixes (missing files, updated documentation, typos) or changes in start or output scripts
    • Medium risk: Uncritical changes in the model core (e.g. moderate modifications in non-default realizations)
    • High risk: Critical changes in model core or default settings (e.g. changing a model default or adjusting a core mechanic in the model)
  • Self-review own code

    • No hard coded numbers and cluster/country/region names.
    • The new code doesn't contain declared but unused parameters or variables.
    • magpie4 R library has been updated accordingly and backwards compatible where necessary.
    • scenario_config.csv has been updated accordingly (important if default.cfg has been updated)
  • Document changes

    • Add changes to CHANGELOG.md
    • Where relevant, put In-code documentation comments
    • Properly address updates in interfaces in the module documentations
    • run goxygen::goxygen() and verify the modified code is properly documented
  • Perform test runs

    • Low risk:
      • Run a compilation check via Rscript start.R --> "compilation check"
    • Medium risk:
      • Run test runs via Rscript start.R --> "test runs"
      • Check logs for errors/warnings
    • High risk:
      • Run test runs via Rscript start.R --> "test runs"
      • Check logs for errors/warnings
      • Default run from the PR target branch for comparison
      • Provide relevant comparison plots (land-use, emissions, food prices, land-use intensity,...)

📉 Performance changes 📈

  • Current develop branch default : 25 mins
  • This PR's default : 25 mins

🚨 Checklist for reviewer 🚨

  • PR is labeled correctly
  • Code changes look reasonable
    • No hard coded numbers and cluster/country/region names.
    • No unnecessary increase in module interfaces
    • model behavior/performance is satisfactory.
  • Changes are properly documented
    • CHANGELOG is updated correctly
    • Updates in interfaces have been properly addressed in the module documentations
    • In-code documentation looks appropriate
  • content review done (at least 1)
  • RSE review done (at least 1)

@flohump flohump added enhancement New feature or request Low risk Low risk labels Nov 6, 2024
@flohump flohump added bugfix and removed Low risk Low risk labels Nov 6, 2024
@flohump flohump marked this pull request as ready for review November 6, 2024 09:29
Copy link
Contributor

@pascal-sauer pascal-sauer left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for fixing these infeasibilities!

v10_balance_positive.scale(j) = 10e-10;
v10_balance_negative.scale(j) = 10e-10;
v10_balance_positive.scale(j,land_from) = 10e-7;
v10_balance_negative.scale(j,land_from) = 10e-7;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you explain what is happening here? Why does it help to switch to 10^-7?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I ran gdx2::calcScaling(gdx) on a couple of different test runs. For most of them the suggestion is to use 10^-7 for scaling.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay this may be more of an understanding question, so the scaling tells gams the order of magnitude of step lengths when iterating towards a local optimum?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

https://www.gams.com/latest/docs/S_CONOPT4.html#CONOPT4_SCALING
internal_model = user_model / scaling_factor
GAMS is also doing automatic scaling but the default for Tol_Scale_Min is 1. Since we have manual scaling factors < 1 I think we should use Tol_Scale_Min = 1.e-10; in the optFile.
https://www.gams.com/latest/docs/S_CONOPT4.html#CONOPTTol_Scale_Min

@@ -23,7 +23,7 @@ else
);
p44_bii_lower_bound(t2,i,biome44)$(p44_bii_lower_bound(t2,i,biome44) >= 1) = 1;
p44_bii_lower_bound(t2,i,biome44)$(m_year(t2) < s44_start_year) = 0;
p44_bii_lower_bound(t2,i,biome44)$(sum(cell(i,j), f44_biome(j,biome44)) = 0) = 0;
p44_bii_lower_bound(t2,i,biome44)$(sum(cell(i,j), f44_biome(j,biome44)) <= 1e-10) = 0;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this change (and other similar changes where 0 is replaced with 1e-10) necessary?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In this specific case for bii_target it was inconsistent with presolve and equations, where it is already 1e-10.
The intention is to fix v44_bii to zero in case the biome share f44_biome is very small. And this case also the lower bound should be zero.
I implemented the same logic for the bii_target_apr24 realisation.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not convinced. Replacing small values with 0 / treating zeros the same as tiny numbers seems like losing control. We don't quite know why this works, but if we replace it with zero it works, and it's almost zero anyway... I have a bad feeling about this practice, we are piling up ultra hard to debug problems in our code this way is my feeling. This is of course more a general comment and not specific to this instance/PR.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand this point. Why should it be an issue to treat tiny numbers like zero?
My understanding is that such tiny numbers might cause issues for the solver.
We also do this in other parts of the model.
For debugging we have off / simple module realizations in many cases.
Do you have another suggestion for addressing this problem?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why should it be an issue to treat tiny numbers like zero?

It's just mathematically wrong, and we don't understand the implications. This lack of control what's going on is what's bugging me, but that is of course inherent to using gams/conopt4, so not something we can change easily :/

We also do this in other parts of the model.

That is not a great argument and just makes me feel more anxious.

Do you have another suggestion for addressing this problem?

Unfortunately no, I cannot follow what is happening really, so I'll approve. Wanted to emphasize my concerns though.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would be interested in the other reviewers perspective on this @tscheypidi @FelicitasBeier @k4rst3ns :) Am I just worrying to much? 🙈

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I totally get your concern. It has been a common practice in the past, though. So maybe something to discuss in the next MAgPIE meeting?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same concern here.

On Florians:

My understanding is that such tiny numbers might cause issues for the solver

... I would not know why theoretically and generally this should be an issue.
We observed that over the years, but do not have a great understanding here.

But I am also see @flohump point: this is not specific to this PR. And I like @FelicitasBeier suggestion bringing this up to the MAgPIE meeting and approving this for now.

@flohump flohump requested a review from pascal-sauer November 7, 2024 10:03

q10_transition_to(j2,land_to) ..
sum(land_from, vm_lu_transitions(j2,land_from,land_to)) =e=
vm_land(j2,land_to);

q10_transition_from(j2,land_from) ..
sum(land_to, vm_lu_transitions(j2,land_from,land_to)) =e=
pcm_land(j2,land_from);
vm_land.l(j2,land_from) + v10_balance_positive(j2,land_from) - v10_balance_negative(j2,land_from);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't know this formulation from equations. Thought that's a presolve/postsolve thing. Are you sure that works?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And in the above two equations, you really want the variable and not the level, right? Just to make sure because the first one (q10_transition_matrix) has changed from pcm_land to vm_land, not vm_land.l

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't know this formulation from equations. Thought that's a presolve/postsolve thing. Are you sure that works?

Yes, it works ;-)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And in the above two equations, you really want the variable and not the level, right? Just to make sure because the first one (q10_transition_matrix) has changed from pcm_land to vm_land, not vm_land.l

I was thinking about the reasons for the infeasible solutions and came to the conclusiong that it might be due to inconsistencies in these 3 equations from the use of vm_land and pcm_land.
Therefore, I switched to vm_land in the first two equations and vm_land.l in third equation.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't fully understand why, but I am happy the infeasibilities are gone, so I'm gonna approve.

Copy link
Member

@k4rst3ns k4rst3ns left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not too much understanding of the changes I believe.
Happy to get feedback on the pcm_land to vm_land.l change before approval (sorry 😅 )


q10_transition_to(j2,land_to) ..
sum(land_from, vm_lu_transitions(j2,land_from,land_to)) =e=
vm_land(j2,land_to);

q10_transition_from(j2,land_from) ..
sum(land_to, vm_lu_transitions(j2,land_from,land_to)) =e=
pcm_land(j2,land_from);
vm_land.l(j2,land_from) + v10_balance_positive(j2,land_from) - v10_balance_negative(j2,land_from);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I understand that vm_land.l should suppose to be the same as pcm_land.
But as you also pointed out pcm_land is actually changed in some resolve functions, so this is not at all identical or not?

I am not 100% sure why we need this balancing, if we could just define a variable with a level that has the correct value of pcm_land?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In this instance it does not make a difference if we use pcm_land or vm_land.l. They should be identical.
The point is that vm_land.l is changed in the presolve based on parameters, which seem to have a different accuracy than the variable.
I tried without the balancing variables. But without them some of the test runs are infeasible.

@@ -23,7 +23,7 @@ else
);
p44_bii_lower_bound(t2,i,biome44)$(p44_bii_lower_bound(t2,i,biome44) >= 1) = 1;
p44_bii_lower_bound(t2,i,biome44)$(m_year(t2) < s44_start_year) = 0;
p44_bii_lower_bound(t2,i,biome44)$(sum(cell(i,j), f44_biome(j,biome44)) = 0) = 0;
p44_bii_lower_bound(t2,i,biome44)$(sum(cell(i,j), f44_biome(j,biome44)) <= 1e-10) = 0;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same concern here.

On Florians:

My understanding is that such tiny numbers might cause issues for the solver

... I would not know why theoretically and generally this should be an issue.
We observed that over the years, but do not have a great understanding here.

But I am also see @flohump point: this is not specific to this PR. And I like @FelicitasBeier suggestion bringing this up to the MAgPIE meeting and approving this for now.

@flohump flohump merged commit 321bea6 into magpiemodel:develop Nov 11, 2024
1 check passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bugfix enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants