Two Stage Stochastic Optimization
for Fixing Energy Reserves

Michel De Lara, Jean-Philippe Chancelier
(last modiﬁcation date: March 18, 2019)

Version pdf de ce document
Version sans bandeaux

1 Problem statement
2 Formulation on a tree with linear costs
3 Formulation on a tree with quadratic convex costs
4 Formulation on a fan with quadratic convex costs
5 Formulation on a fan with linear costs
A Additional code for “Formulation on a tree with linear costs”
B Additional code for “Formulation on a tree with quadratic convex costs”

1 Problem statement

We formalize the problem of ﬁxing energy reserves in a day-ahead market as a two stage stochastic optimization problem. A decision has to be made at night of day J — which quantity of the cheapest energy production units (reserve) to be mobilized — to meet a demand that will materialize at morning of day J + 1. Excess reserves are penalized; demand unsatisﬁed by reserves have to be covered by costly extra units. Hence, there is a trade-oﬀ to be assessed by optimization.

1.1 Stages

There are two stages, represented by the letter t (for time):

t = 0 corresponds to night of day J;
t = 1 corresponds to morning of day J + 1.

1.2 Probabilistic model

We suppose that the demand, materialized on the morning of day J + 1, can take a ﬁnite number S of possible values D_s, where s denotes a scenario in the ﬁnite set 𝕊 (S=card(𝕊)).

We denote π_s the probability of scenario s, with

∑ ∀s ∈ 𝕊, πs > 0, πs = 1. s∈𝕊

(1)

Notice that we do not consider scenarios with zero probability.

1.3 Decision variables

The decision variables are the scalar Q₀ and the ﬁnite sequence (Q1,s) _s∈𝕊 of scalars, as follows:

at stage t = 0, the energy reserve is Q₀;
at stage t = 1, a scenario s materializes and the demand D_s is observed, so that one decides of the recourse quantity Q_1,s.

The decision variables can be considered as indexed by a tree with one root (corresponding to the index 0) and as many leafs as scenarios in 𝕊 (each leaf corresponding to the index 1,s): Q₀ is attached to the root of the tree, and each Q_1,s is attached to a leaf corresponding to s.

1.4 Optimization problem formulation

The balance equation between supply and demand is

Q0 + Q1,s = Ds, ∀s ∈ 𝕊.

(2)

The energies mobilized at stages t = 0 and t = 1 display diﬀerent features:

at stage t = 0, the energy production has maximal capacity Q₀^♯, and producing Q₀ costs 𝒞₀(Q₀);
at stage t = 1, the energy production is supposed to be unbounded, and producing Q₁ costs 𝒞₁(Q₁).

We consider the stochastic optimization problem

Here, we look for energy reserve Q₀ and recourse energy Q_1,s so that the balance equation (3d) is satisﬁed (at stage t = 1) at minimum expected cost in (3a). By weighing each scenario s with its probability π_s, the optimal solution (Q⋆,(Q ⋆ ) )
0 1,s s∈𝕊 performs a compromise between scenarios.

2 Formulation on a tree with linear costs

Here, we suppose that the costs are linear:

𝒞0 (Q0 ) = l0Q0, 𝒞1(Q1 ) = l1Q1.

(4)

Therefore, the stochastic optimization problem (3) now becomes

This optimization problem (5) is linear. When the number S of scenarios is not too large, we can use linear solvers.

Question 1 We consider the case when 𝕊 = {L,M,H} has S = 3 scenarios (low, medium, high). We want to transform the linear optimization problem (5) under a form adapted to a linear solver.

a): [1+1+1] Expand the criterion (5a). Expand the inequalities (5b)–(5c) into an array of ﬁve scalar equations (one inequality per equation), and the equalities (5d) into an array of three scalar equations (one equality per equation).
b): [1+1+1+1+1] Let x be the column vector x = (Q₀,Q_1L,Q_1M,Q_1H). Propose matrices A_e, A_i and column vectors b_e, b_i and c such that the linear optimization problem (5) can be written under the form

Question 2 We are going to numerically solve the linear optimization problem (5).

a): [2] Interpret the code below. What is the macro linpro doing? What is lopt? Copy the code into a ﬁle named tp_q1.sce.
b): [1+2] Solve a numerical version of problem (5) with S = 3 scenarios and the parameters in the code below, by executing the ﬁle tp_q1.sce. What is the optimal value Q₀^⋆ of the reserve? What is the optimal value Q_1,L^⋆? Can you explain why these values are optimal? (there is an economic explanation based on relative costs; there is a mathematical explanation based on the properties of the solutions of a linear program).
c): [3] Then, increase and decrease the value of the unitary cost l₁ (especially above and below l₀). Show the numerical results that you obtain. What happens to the optimal values Q₀^⋆? Explain why (make the connection with the properties of the solutions of a linear program).

  // Formulation on a tree with linear costs.
  // Numerical resolution by linear programming

  // Constant initialization
  S=3;// Number of random scenarios
  q0m=30;// max capacity for q0

  // Demand
  if S==3 then
    D=[15;20;50];
    Pr=[0.2;0.6;0.2];// Probabilities of Demand
  else
    D=grand(S,1,'uin',5,50);
    Pr=grand(S,1,'unf',0,1);
    Pr=Pr ./sum(Pr);// Probabilities of Demand
  end

  // Constants used in the cost function
  ll0=2;ll1=5;

  // a revoir pour passer a des contraintes egalité et utiliser lb et ub
  c=[ll0;ll1 .*Pr];// cost coefficients

  // inequality constraints (bounds on production)
  Ai=[-eye(S+1,S+1);eye(1,S+1)];
  bi=[zeros(S+1,1);q0m];

  // equality constraints, i.e. production equals demand
  Ae=[ones(S,1),eye(S,S)];
  be=[D];

  // solving by linear
  // xopt should be [  15,  0, 5, 35 ] when S=3

  // scicoslab version with linpro
  A=[Ae;Ai];b=[be;bi];
  [xopt,lopt,fopt]=linpro(c,A,b,[],[],size(be,'*'))

Question 3 We are going study the impact of the number S of scenarios on the numerical resolution of the linear optimization problem (5).

a): [2] Take S = 100. What is the optimal value Q₀^⋆ of the reserve? Identify the scenario s with the lowest demand. What is the optimal value Q_1,s^⋆? Explain.
b): [1] For what value of n in S = 10ⁿ can you no longer solve numerically?

3 Formulation on a tree with quadratic convex costs

Here, we suppose that the costs are quadratic and convex:

𝒞0(Q0) = 1-K0Q2 + l0Q0, K0 > 0, 𝒞1(Q1 ) = 1K1Q2 + l1Q1, K1 > 0. 2 0 2 1

(7)

The optimization problem (3) is quadratic convex:

When the number S of scenarios is not too large, we can use quadratic solvers.

Question 4 We are going to numerically solve the quadratic convex optimization problem (8).

a)

[1] Interpret the code below. What is the macro quapro doing? Copy the code into a ﬁle named tp_q2.sce.

b)

[1] Solve a numerical version with S = 3 scenarios and the parameters in the code below, by executing the ﬁles tp_q1.sce and tp_q2.sce. What is the optimal value Q₀^⋆ of the reserve? What is the optimal value Q_1,L^⋆?

c)

[1+1+2] Check numerically that $(Q ⋆0,Q ⋆1,L,Q ⋆1,H ,Q ⋆1,M )$ is an inner solution to the optimization problem (8a), that is, check numerically that the inequalities (8b) and (8c) are strict. What is the diﬀerence with the optimal solution of Question 2? Discuss the diﬀerence (make the connection with the properties of the solutions of a linear program).

d)

[2+1] Compute the derivatives of the cost functions in (7). Check numerically (giving the details of computation) that the optimal solution $(Q ⋆0,Q ⋆1,L,Q ⋆1,M ,Q ⋆1,H)$ satisﬁes the following relation between marginal costs:

𝒞′0(Q ⋆0) = πL𝒞′1(Q ⋆1,L) + πM 𝒞′1(Q ⋆1,M ) + πH 𝒞 ′1(Q ⋆1,H ).

(9)

  // Formulation on a tree with quadratic costs
  // Numerical resolution by brute force quadratic programming

  // just to get common data
  exec('tp_q1.sce');

  // the problem is now quadratic; we use a quadratic solver
  A=[Ae;Ai];b=[be;bi];
  KK0=10;KK1=1;
  KK=diag([KK0,KK1*Pr' .*ones(1,S)]);

  // scicoslab version with quapro
  [xopt,lopt,fopt]=quapro(KK,c,A,b,[],[],size(be,'*'))

Question 5 We are going to numerically solve the quadratic convex optimization problem (8) after changing the relative values of the (quadratic) parameters K₀ and K₁. For the parameters K₀, l₀, l₁, we take the same values as those in Question 2 (hence l₁ > l₀) and in Question 4.

a): [1] Take K₁ > K₀ with K₁ ≈ K₀. What are the optimal value Q₀^⋆ of the reserve and the optimal value Q_1,L^⋆?
b): [1] Take K₁ > K₀ with K₁ >> K₀. What are the optimal value Q₀^⋆ of the reserve and the optimal value Q_1,L^⋆?
c): [2] Discuss.

Question 6 We are going study the impact of the number S of scenarios on the numerical resolution of the quadratic convex optimization problem (8).

a): [1] Take S = 100. Solve a numerical version of problem (8) with the parameters K₀, l₀, K₁, l₁ as in the code of Question 4. What is the optimal value Q₀^⋆ of the reserve? Identify the scenario s with the lowest demand. What is the optimal value Q_1,s^⋆?
b): [1] For what value of n in S = 10ⁿ can you no longer solve numerically?

Question 7 This theoretical question may be ignored
(by those who want to focus on numerical results)
For this question, we temporarily ignore the inequalities (8b) and (8c) in (8). Therefore, we consider the optimization problem (8a) with equality constraints (8d), that is:

a)

[2] Compute the Hessian matrix of the criterion

[ ] S ∑ 1- 2 1- 2 J0 : (Q0, (Q1,s)s∈𝕊) ∈ ℝ × ℝ ↦→ πs 2K0Q 0 + l0Q0 + 2 K1Q 1,s + l1Q1,s . s∈𝕊

(11)

What are the dimensions of the Hessian matrix?

b)

[3] Why does optimization problem (10) have a solution? (Beware of the domain)

c)

[2] Why is the solution unique?

d)

[1+2] Why are the equality constraints (10b) qualiﬁed? Why does an optimal solution $⋆ ⋆ (Q 0,(Q 1,s)s∈𝕊)$ of (10) satisfy the Karush-Kuhn-Tucker (KKT) conditions (ﬁrst-order optimality conditions)?

e)

[2] Why is a solution of the KKT conditions an optimal solution of (10)?

f)

[1] Write the Lagrangian $ℒ0 (Q0, (Q1,s)s∈𝕊,(μs)s∈𝕊)$ associated with problem (10).

g)

[2] Deduce the KKT conditions. Show that there exist (μ_s^⋆,s ∈ 𝕊) such that

′ ⋆ ∑ ⋆ ′ ⋆ ⋆ 𝒞0(Q 0) − μs = 0 and πs𝒞 1(Q 1,s) − μs = 0, ∀s ∈ 𝕊. s∈ 𝕊

(12)

h)

[2] Deduce that — when $(Q ⋆0,(Q ⋆1,s)s∈ 𝕊)$ is an inner optimal solution to problem (8a) — we have the following relation between marginal costs:

∑ 𝒞′(Q ⋆) = π 𝒞 ′(Q ⋆ ). 0 0 s 1 1,s s∈ 𝕊

(13)

Give an economic interpretation of this equality.

4 Formulation on a fan with quadratic convex costs

When the number S of scenarios is too large, Problem (5) — be it linear or convex — cannot be solved by direct methods.

4.1 Dualization of non-anticipativity constraints

To bypass this problem, we use a “trick” consisting in introducing new decision variables (Q0,s) _s∈𝕊, instead of the single decision variable Q₀, and we write

∑ Q0,s = πs′Q0,s′, ∀s ∈ 𝕊. s′∈ 𝕊

(14)

These equalities are called the non-anticipativity constraints. Indeed, the equations (14) express that

Q = Q ′, ∀(s,s′) ∈ 𝕊2, 0,s 0,s

(15)

that is, the decision at stage t = 0 does not depend on the scenario s, hence cannot anticipate the future. Later, we will treat the constraints (14) by duality.

Therefore, the stochastic optimization problem (3) now becomes

By the assumption (1) that there are no scenarios with zero probability (π_s > 0), we replace each equality in the last equation by the equivalent one

∑ πsQ0,s = πs πs′Q0,s′, ∀s ∈ 𝕊. s′∈𝕊

(16e)

We attach, to each equality above a multiplier λ_0,s. We put

Q = ((Q0,s)s∈𝕊,(Q1,s)s∈𝕊), λ = (λ0,s) . s∈𝕊

(17)

The corresponding Lagrangian is

Question 8 This theoretical question may be ignored
(by those who want to focus on numerical results)
When the costs are quadratic and convex as in (7), show that

a)

[3] the criterion

[ ] S S ∑ 1 2 1 2 J : Q ∈ ℝ × ℝ ↦→ πs 2K0Q 0,s + l0Q0,s + 2-K1Q 1,s + l1Q1,s s∈𝕊

(19)

in (16a) is a-strongly convex in $((Q0,s)s∈𝕊,(Q1,s)s∈𝕊)$ , and provide a possible a > 0,

b)

[2] the domain deﬁned by the constraints in the optimization problem (16) is closed,

c)

[2] the domain deﬁned by the constraints in the optimization problem (16) is convex,

d)

[2+1] the optimization problem (16) has a solution, and it is unique, denoted by Q^⋆,

e)

[1] there exists a multiplier λ^⋆ such that (Q^⋆,λ^⋆) is a saddle point of the Lagrangian ℒ in (18).

4.2 Uzawa algorithm

We consider the following optimization problem

under the assumptions that

the set 𝕌^ad is a closed convex subset of a Euclidian space ℝ^N,
the criterion J : ℝ^N → ℝ is an a-strongly convex (a > 0) and diﬀerentiable function,
the constraint mapping 𝜃 : ℝ^N → ℝ^M is aﬃne, κ-Lipschitz (κ > 0),
the Lagrangian ℒ(u,λ) = J(u) + ⟨λ, Θ(u)⟩ admits a saddle-point over 𝕌^ad × ℝ^M.

Then the following algorithm — called dual gradient algorithm, or Uzawa algorithm — converges toward the optimal solution of Problem (20), when 0 < ρ < 2a∕κ².

Data: Initial multiplier λ⁽⁰⁾, step ρ Result: optimal decision and multiplier; repeat u^(k) = arg min _u∈𝕌^adℒ(u,λ^(k)) ;
λ^(k+1) = λ^(k) + ρΘ(u^(k)) ; until Θ(u^(k)) = 0; Algorithm .1: Dual Gradient Algorithm

4.3 Numerical resolution by Uzawa algorithm (quadratic convex case)

When the costs are quadratic and convex as in (7), the optimization problem (16) becomes:

Question 9 This theoretical question may be ignored
(by those who want to focus on numerical results)
When the costs are quadratic and convex as in (7), identify in the optimization problem (21) the corresponding elements in the Uzawa algorithm .1:

a): [1] decision variable u,
b): [1] decision set ℝ^N (for what N?),
c): [2] aﬃne constraints mapping Θ : ℝ^N → ℝ^M, corresponding to the constraints (21e) (why is it κ-Lipschitz, and for which κ?).
d): [3] Explain why the Uzawa algorithm converges towards the optimal solution of Problem (21).

Question 10 We are going to numerically solve the quadratic convex optimization problem (21) by Uzawa algorithm.

a): [3] Detail what the code below is doing ; explain how the code implements the Uzawa algorithm. Why do we use the macro quapro? What is the dual function? Copy the code into a ﬁle named tp_q3.sce.
b): [2] Solve a numerical version with S = 3 scenarios. What is the solution $(Q ⋆,(Q ⋆ ) ) 0 1,ss∈𝕊$ given by the algorithm? Do you have that $D = Q ⋆+ Q⋆ 1,s 0 1,s$ , for all $s ∈ 𝕊$ ?
c): [1] Then, try with S = 100. What is the value Q₀^⋆ of the reserve given by the algorithm? Identify the scenario s with the lowest demand. What is the value Q_1,s^⋆ given by the algorithm?
d): [1] For what value of n in S = 10ⁿ can you no longer solve numerically?

  // formulation on a fan
  // Constant initialization

  // Constant initialization
  S=3;// Number of random scenarios
  q0m=30;// max capacity for q0

  // Demand
  if S==3 then
    D=[15;20;50];
    Pr=[0.2;0.6;0.2];// Probabilities of Demand
  else
    D=grand(S,1,'uin',5,50);
    Pr=grand(S,1,'unf',0,1);
    Pr=Pr ./sum(Pr);// Probabilities of Demand
  end

  // Constants used in the cost function
  ll0=2;ll1=5;
  KK0=10;KK1=1;

  // Uzawa iterations when the dualized constraints are the S constraints
  //  Q0(ss) = sum(Pr.*Q0);

  Q0=zeros(S,1);
  Q1=zeros(S,1);
  f0=zeros(S,1);

  rho=5;

  lambda=zeros(S,1);
  // iterations of the Uzawa algorithm
  for it=0:30 do
    // decomposed minimizations (loop over scenarios ss)
    for ss=1:S do
      // inequality constraints (bounds)
      // bounds on (Q0s,Q1s)
      ll=[ll0*Pr(ss);ll1*Pr(ss)];// cost coefficients
      Ai=[-eye(2,2);eye(1,2)];
      bi=[zeros(2,1);q0m];
      // equality constraints i.e production equals demand
      Ae=[1,1];
      be=[D(ss)];
      A=[Ae;Ai];b=[be;bi];
      KK=Pr(ss)*diag([KK0,KK1]);
      //
      cc=ll+Pr(ss)*[lambda(ss)-sum(Pr .*lambda);0];
      [xopt,lopt,fopt]=quapro(KK,cc,A,b,[],[],size(be,'*'));
      Q0(ss)=xopt(1);
      Q1(ss)=xopt(2);
      f0(ss)=fopt;
      printf("Dual function %f\n",sum(f0));
    end
    Q0bar=sum(Pr .*Q0);
    lambda=lambda+rho*(Pr .*(Q0-Q0bar));
  end

5 Formulation on a fan with linear costs

Here, we suppose that the costs are linear, as in (4).

5.1 Diﬃculties in applying Uzawa algorithm with linear costs

The optimization problem (16) becomes:

Question 11 We are going to numerically solve the linear optimization problem (22).

a): [2] Detail what the code below is doing. Why do we use the macro linpro? Copy the code into a ﬁle named tp_q4.sce.
b): [3] Solve a numerical version with S = 3 scenarios. What do you observe regarding convergence of the Uzawa algorithm? Can you explain why?

  // formulation on a fan with linear cost
  // Uzawa does not work

  S=3;// Number of random scenarios
  q0m=30;// max capacity for q0
  D=[15;20;50];
  Pr=[0.2;0.6;0.2];// Probabilities of Demand
  rho=0.1;
  // Constants used in the cost function
  ll0=2;ll1=5;

  // Uzawa iterations when the dualized constraints are the S constraints
  //  Q0(ss) = sum(Pr.*Q0);

  Q0=zeros(S,1);
  Q1=zeros(S,1);
  f0=zeros(S,1);

  lambda=zeros(S,1);
  for it=0:30 do
    // decomposed minimizations
    for ss=1:S do
      // inequality constraints (bounds)
      // bounds on (Q0s,Q1s)
      ll=[ll0*Pr(ss);ll1*Pr(ss)];// cost coefficients
      Ai=[-eye(2,2);eye(1,2)];
      bi=[zeros(2,1);q0m];
      // equality constraints i.e production equals demand
      Ae=[1,1];
      be=[D(ss)];
      A=[Ae;Ai];b=[be;bi];
      //
      cc=ll+Pr(ss)*[lambda(ss)-sum(Pr .*lambda);0];
      [xopt,lopt,fopt]=linpro(cc,A,b,[],[],size(be,'*'));
      Q0(ss)=xopt(1);
      Q1(ss)=xopt(2);
      f0(ss)=fopt;
      printf("Dual function %f\n",sum(f0));
    end
    Q0bar=sum(Pr .*Q0);
    lambda=lambda+rho*(Pr .*(Q0-Q0bar));
  end

5.2 Dualization of non-anticipativity constraints

To bypass the diﬃculty in applying Uzawa algorithm with linear costs, we use a “trick” consisting in introducing new decision variables Q₀ and (Q )
0,s _s∈𝕊, instead of the single decision variable Q₀, and we write

--- Q0,s = Q0, ∀s ∈ 𝕊.

(23)

These equalities are another form of the non-anticipativity constraints. Indeed, the equations (23) express that the decision at stage t = 0 cannot anticipate the future, hence cannot depend on the scenario s. We will treat these constraints by duality.

Therefore, the stochastic optimization problem (3) now becomes

By the assumption (1) that there are no scenarios with zero probability (π_s > 0), we replace each equality in (24e) by the equivalent one

--- πsQ0,s = πsQ0, ∀s ∈ 𝕊.

(25)

We attach, to each equality above a multiplier λ_0,s. We put

Q0 = (Q0,s) , Q1 = (Q1,s) , λ = (λ0,s) . s∈ 𝕊 s∈𝕊 s∈𝕊

(26)

5.3 Augmented Lagrangian and obstacles to decomposition

5.4 Progressive Hedging algorithm (quadratic solver)

Question 12 We are going to numerically solve the linear optimization problem (24) by the Progressive Hedging algorithm.

a): [3] Detail what the code below is doing. Why do we use the macro quapro? Explain the two roles of the new parameter rr. Copy the code into a ﬁle named tp_q5.sce.
b): [2] Solve a numerical version with S = 3 scenarios. What do you observe regarding convergence of the Uzawa algorithm? Can you explain why?

  // formulation on a fan
  // with linear cost and augmented lagrangian

  S=3;// Number of random scenarios
  q0m=30;// max capacity for q0
  D=[15;20;50];
  Pr=[0.2;0.6;0.2];// Probabilities of Demand
  // Constants used in the cost function
  ll0=2;ll1=5;
  // Constant used both as a quadratic term and as a gradient step
  rr=0.1;

  // Uzawa iterations when the dualized constraints are the S constraints
  //   Q0s = sum(Pr.*Q0)

  Q0=zeros(S,1);
  Q1=zeros(S,1);
  f0=zeros(S,1);

  lambda=zeros(S,1);
  Q0bar=0;
  for it=0:300 do
    // decomposed minimizations
    // we alternate minimization
    for ss=1:S do
      // inequality constraints (bounds)
      // bounds on (Q0s,Q1s)
      ll=[ll0*Pr(ss);ll1*Pr(ss)];// cost coefficients
      Ai=[-eye(2,2);eye(1,2)];
      bi=[zeros(2,1);q0m];
      // equality constraints i.e production equals demand
      Ae=[1,1];
      be=[D(ss)];
      A=[Ae;Ai];b=[be;bi];
      KK=rr*Pr(ss)*diag([1,0]);
      cc=ll+[Pr(ss)*lambda(ss);0];
      cc=cc+[-rr*Q0bar*Pr(ss)];
      [xopt,lopt,fopt]=quapro(KK,cc,A,b,[],[],size(be,'*'));
      Q0(ss)=xopt(1);
      Q1(ss)=xopt(2);
      f0(ss)=fopt;
    end
    // updates of Q0bar
    Q0bar=sum(Pr .*Q0);
    lambda=lambda+rr*(Pr .*(Q0-Q0bar));
  end

  // solution (Q0bar,Q1)
  Q1=max(0,(D-Q0bar));
  // to be compared with tp_q1

5.5 Progressive Hedging algorithm (linear solver)

A Additional code for “Formulation on a tree with linear costs”

  // Formulation on a tree with linear costs.
  // Numerical resolution by linear programming
  // using nsp linprog (glpk).

  // Constant initialization
  S=3;// Number of random scenarios
  q0m=30;// max capacity for q0

  // Demand
  if S==3 then
    D=[15;20;50];
    Pr=[0.2;0.6;0.2];// Probabilities of Demand
  else
    D=grand(S,1,'uin',5,50);
    Pr=grand(S,1,'unf',0,1);
    Pr=Pr ./sum(Pr);// Probabilities of Demand
  end

  // Constants used in the cost function
  ll0=2;ll1=5;

  // a revoir pour passer a des contraintes egalité et utiliser lb et ub
  c=[ll0;ll1 .*Pr];// cost coefficients

  // inequality constraints (bounds on production)
  Ai=[-eye(S+1,S+1);eye(1,S+1)];
  bi=[zeros(S+1,1);q0m];

  // equality constraints, i.e. production equals demand
  Ae=[ones(S,1),eye(S,S)];
  be=[D];

  // solving by linear
  // xopt should be [  15,  0, 5, 35 ] when S=3

  if exists('%nsp') then
    // in nsp linprog is glpk
    [xopt,fopt,flag]=linprog(c,Ai,bi,Ae,be);
    // we can also use quapro if available with a 0 quadratic term
    if exists('quapro','callable') then
      A=[Ae;Ai];b=[be;bi];
      [xopt_q,lagr_q,fopt_q]=quapro(zeros(S+1,S+1),c,A,b,[],[],size(be,'*'))
    end
  else
    // scicoslab version with linpro
    A=[Ae;Ai];b=[be;bi];
    [xopt,lopt,fopt]=linpro(c,A,b,[],[],size(be,'*'))
  end

B Additional code for “Formulation on a tree with quadratic convex costs”

  // Formulation on a tree with quadratic costs
  // Numerical resolution by brute force quadratic programming
  // using nsp quapro or optim

  if exists('%nsp') then
    load_toolbox('quapro');
  end

  // just to get common data
  exec('tp_q1.sce');

  // the problem is now quadratic; we use a quadratic solver if available
  A=[Ae;Ai];b=[be;bi];
  KK0=10;KK1=1;
  KK=diag([KK0,KK1*Pr' .*ones(1,S)]);

  if exists('%nsp') then
    // we can also use quapro if available with a 0 quadratic term
    if exists('quapro','callable') then
      [xopt_q,lagr_q,fopt_q]=quapro(KK,c,A,b,[],[],size(be,'*'))
    end
  else
    // scicoslab version with quapro
    [xopt,lopt,fopt]=quapro(KK,c,A,b,[],[],size(be,'*'))
  end

  // BEWARE: using a global optim if quapro is not available
  // we eliminate the equality constraint to use optim on q0

  function [f,g,ind]=costf(q0,ind)
    q1=max(0,(D-q0));
    // expression of the cost
    f=sum(Pr .*(ll0*q0+ll1*q1))+(1/2)*[q0;q1]'*KK*[q0;q1]
    // expression of the cost derivative
    g=sum(Pr .*(ll0-ll1*(D-q0 >= 0)))+KK(1,1)*q0+-sum((D-q0 >= 0) .*(KK(2:$,2:$)*q1))
  endfunction

  if exists('%nsp') then
    Q0=optim(costf,0,xinf = 0,xsup = q0m)
  else
    [fo,Q0]=optim(costf,'b',0,q0m,0)
  end

  // solution (Q0,Q1)
  Q1=max(0,(D-Q0))

Formations

Départements

Pratique

Recherche

Management of Energy

Management of Renewable Resources and Biodiversity

Management of Exhaustible Resources

Models in Life Science

Miscellaneous

Two Stage Stochastic Optimization
for Fixing Energy Reserves

1 Problem statement

1.1 Stages

1.2 Probabilistic model

1.3 Decision variables

1.4 Optimization problem formulation

2 Formulation on a tree with linear costs

3 Formulation on a tree with quadratic convex costs

4 Formulation on a fan with quadratic convex costs

4.1 Dualization of non-anticipativity constraints

4.2 Uzawa algorithm

4.3 Numerical resolution by Uzawa algorithm (quadratic convex case)

5 Formulation on a fan with linear costs

5.1 Diﬃculties in applying Uzawa algorithm with linear costs

5.2 Dualization of non-anticipativity constraints

5.3 Augmented Lagrangian and obstacles to decomposition

5.4 Progressive Hedging algorithm (quadratic solver)

5.5 Progressive Hedging algorithm (linear solver)

A Additional code for “Formulation on a tree with linear costs”

B Additional code for “Formulation on a tree with quadratic convex costs”

L'École des Ponts ParisTech est membre fondateur de	L'École des Ponts ParisTech est certifiée

Formations

Départements

Pratique

Recherche

Management of Energy

Management of Renewable Resources and Biodiversity

Management of Exhaustible Resources

Models in Life Science

Miscellaneous

Two Stage Stochastic Optimization for Fixing Energy Reserves

1 Problem statement

1.1 Stages

1.2 Probabilistic model

1.3 Decision variables

1.4 Optimization problem formulation

2 Formulation on a tree with linear costs

3 Formulation on a tree with quadratic convex costs

4 Formulation on a fan with quadratic convex costs

4.1 Dualization of non-anticipativity constraints

4.2 Uzawa algorithm

4.3 Numerical resolution by Uzawa algorithm (quadratic convex case)

5 Formulation on a fan with linear costs

5.1 Diﬃculties in applying Uzawa algorithm with linear costs

5.2 Dualization of non-anticipativity constraints

5.3 Augmented Lagrangian and obstacles to decomposition

5.4 Progressive Hedging algorithm (quadratic solver)

5.5 Progressive Hedging algorithm (linear solver)

A Additional code for “Formulation on a tree with linear costs”

B Additional code for “Formulation on a tree with quadratic convex costs”

Two Stage Stochastic Optimization
for Fixing Energy Reserves