Bidding System Design for Multiperiod Electricity Markets:
Pricing of Stored Energy Shiftability
Takayuki Ishizaki1 , Masakazu Koike2 , Nobuyuki Yamaguchi3 , and Jun-ichi Imura1
Abstract— In this paper, we design a bidding system fora multiperiod electricity market in which market players
participate with power generators and energy storage resources.
For the bidding system design, we first develop a sequential procedure to determine a separate multidimensional bid
function, i.e., an ensemble of period-specific bid functions,which enables to regard the multiperiod electricity market
as an ensemble of conventional period-specific electricity markets. This sequential determination also enables to construct
a distributed approximate scheme for multiperiod electricity
market clearing. Then, based on a basis transformation similar
to the Fourier transformation, we propose a bidding systemwith explicit consideration of the pricing of energy shiftability.
It is shown that, in the situation where the optimal priceprofile levels off due to high penetration of energy storage, thedistributed approximate scheme in the Fourier-like basis canattain the optimal market clearing with the minimal social cost.
In addition, we numerically investigate the resultant deadweight
loss, i.e., an increase of social costs caused by approximation,
varying the levels of energy storage penetration.

I. I NTRODUCTION
The development of a smart grid has been recognized as
one of key issues in addressing environmental and social concerns; see [1], [2] for pedagogical overviews. In particular,
towards effective integration of dispatchable and renewable
power generation, the potential of energy storage has been
attracting international attention in smart grid community.
Along this trend, developing a multiperiod electricity market, as opposed to a conventional period-specific electricity
market, is crucial for making use of the power shiftability ofenergy storage resources and flexible loads. The significanceand difficulty of developing such an electricity market has
been discussed in the literature; see, e.g., Section 4.3 of [1].
With this background, we have formulated in [3] a multiperiod electricity market as a generalization of a conventional
period-specific electricity market. In this formulation, eachaggregator aims at making the highest profit by transacting
the prosumption (production and consumption) of energy
amounts at multiple timeslots, called a prosumption profile. The prosumption profile is represented as a vector of
prosumption energy amounts. Each aggregator generates a
1 Department of Mechanical and Environmental Informatics, Graduate
School of Information Science and Engineering, Tokyo Institute of Technology; 2-12-1, Meguro, Tokyo, Japan.
2 Department of Marine Electoronics and Mechanical Engineering, Tokyo
University of Marine Science and Technology; 4-5-7, Kounan, Minato,
Tokyo, 108-0075, Japan.
3 Tokyo University of Science; 6-3-1, Niijuku, Katsushika, Tokyo 1258585, Japan.
This research was supported by CREST, JST.
{ishizaki,imura}@mei.titech.ac.jp
mkoike0@kaiyodai.ac.jp

prosumption profile based on available energy resources,
such as dispatchable power generation and energy storage.
Devising a clearing scheme for the multiperiod electricity
market is not straightforward due to the multidimensionality(time dependence) of prosumption profiles; see Section IIIA for the details. For this issue, we have also proposed in[3] a distributed market clearing scheme implemented as
an indirect communication among aggregators through an
independent system operator (ISO). This is developed basedon the premise that the optimal clearing price profile exactly
levels off. Even though we give a theoretical clarificationthat such exact price levelling off is expected when thepenetration of energy storage is sufficiently high, devisinga market clearing scheme without assuming the exact pricelevelling off is left as an open question there.
As generalization of the previous work, this paper develops
a more general bidding system for the multiperiod electricity
market that can attain prosumption balancing without assuming the exact price levelling off. To this end, we discuss thefollowing items:– We propose a sequential procedure to determine an
ensemble of conventional period-specific bid functions
for the multiperiod electricity market.
– Based on a basis transformation similar to the Fouriertransformation, we propose a bidding system that takes
explicit account of the pricing of energy shiftability.
Section IV in this paper is devoted to the first item.
The proposed sequential procedure can be seen as an approximation method to reduce an additively indecomposable
cost function to an additively decomposable function. In
general, the cost function of energy storage is not additively decomposable because it is temporally dependent.
The additively decomposable approximant leads to an ensemble of period-specific bid functions, which enables toregard the multiperiod electricity market as an ensemble ofconventional period-specific electricity markets. This alsoenables to construct a bidding system implemented as a
distributed approximate scheme for market clearing. It should
be noted that this bidding system necessarily complies withthe constraint of prosumption profile balance even though it
may cause a degree of deadweight losses.
Section V is devoted to the second item. The transformed
basis, called the multiresolved basis, represents an orderedbasis compatible with temporal resolution (frequency) of
prosumption and price profiles. The corresponding bidding
system in this basis is regarded as a situation where a market
for the total prosumption amount on the day of interest opens
firstly, a market for shifting prosumption amounts between

the morning timeslots and the afternoon timeslots opens
secondly, and so force. This pricing of energy shiftability isnovel compared with the existing works [4]–[7]. In addition,
it is clarified that market clearing with the minimal social
cost, i.e., no deadweight loss, can be attained under highpenetration of energy storage.
The remainder of this paper is organized as follows. We
first overview the formulation of the multiperiod electricity
market in Section II. Then, in Section III, we formulate
a bidding system design problem. The main contributions
are provided in Sections IV and V, both of which include
illustrative numerical examples. Finally, concluding remarks
are provided in Section VI.
Notation: We denote the set of real values by R, the ith
column of the identity matrix by ei , the n-dimensional allones vector by 1n , and the direct product of S1 , . . . , Sn by∏
S1 × · · · × Sn = i∈{1,...,n} Si .
A function F : Rn → R is said to be convex if()
F (1 − λ)x + λx′ ≤ (1 − λ)F (x) + λF (x′ )for every pair of x and x′ in X . An interval domain ofx ∈ Rn is denoted by [x] ⊂ Rn .
II. M ULTIPERIOD E LECTRICITY M ARKETS
A. Aggregator ModelFirst, we overview the multiperiod electricity market [3]in which several aggregators participate. We provide a mathematical model of one aggregator who provides the prosumption (production and consumption) of energy amounts
on the day of interest. Let T := {1, 2, . . . , n} denote theset of timeslots on the day. The prosumption energy of anaggregator at the tth timeslot can be described ast∈T

Furthermore, we denote the set of realizable prosumption
profiles as{}
X := x ∈ Rn : F(x) ̸= ∅ ,(4)which is convex if G and D are convex. The cost functionsof dispatchable power generation and battery charge anddischarge are denoted as
G : G → R,

(1)

for all λ ∈ (0, 1) and for every pair of x and x′ in the domainsuch that the value of F is finite. In particular, F is said tobe strictly convex if (1) holds with the strict inequality unlessx = x′ . A set-valued function f : Rn → Rn is said to be
monotone increasing over X if()Ty − y ′ (x − x′ ) ≥ 0, ∀y ∈ f (x), y ′ ∈ f (x′ )xt = gt − lt + η out δtout − η1in δtin ,charge and discharge power profiles δ in and δ out are decisionvariables. To realize a desired prosumption profile x, theaggregator determines g and δ := (δ in , δ out ) as complying
with the constraints of g ∈ G and δ ∈ D, which represents thebounds for dispatchable power generation and the limitationof inverter and battery capacities, respectively.
With respect to each prosumption profile x, we denote thefeasible subspace of the dispatchable power generation andthe battery charge and discharge profiles as{}
F(x) := (g, δ) ∈ G × D : (2) holds .
(3)

(2)

where xt ∈ R denotes the resultant prosumption energy
to the grid, gt ∈ R+ denotes the power generation of
dispatchable generators, lt ∈ R+ denotes the load, andδtin ∈ R+ and δtout ∈ R+ denote the battery charge anddischarge power. The positive constants η in and η out denotethe charge and discharge efficiency, respectively, each ofwhich takes a value in (0, 1]. Note that the sign of xt is
positive for outflow direction to the grid.
In the following, we denote the stacked vector of a symbol,
called a profile, by that without the subscript t. For example,
the prosumption profile x = (xt )t∈T represents the sequenceof prosumption amounts on the day. In the following, the load
profile l is supposed to be a given vector. On the other hand,the dispatchable power generation profile g, and the battery

D : D → R.

(5)

As shown in [3], if G and D are both convex over G and D,
respectively, then the prosumption cost function defined by{}
F (x) := min
G(g) + D(δ)(6)(g,δ)∈F (x)is convex over X . The value of F (x) represents the minimum
cost to realize a prosumption profile x. Note that the closed
form of F cannot be written down in general. Furthermore, F
is generally not strictly convex even though it is necessarily
convex if G and D are convex. This is due to the fact that thebattery charge and discharge cost function D is not strictly
convex because energy storage has the shiftability of energy
amounts among timeslots.
B. Derivation of Multidimensional Bid Functions
In this subsection, we discuss a bid function compatible
with the multiperiod electricity market. Let λ ∈ [λ] denote ann-dimensional price profile for prosumption profile transaction, which corresponds to the sequence of prices on the dayof interest. The profit function in transacting a prosumption
profile x is defined as
J(x; λ) := λT x − F (x),(7)

where λT x represents the total income by transacting the
prosumption profile. On the basis of this profit function, theaggregator can determine a set-valued function given as{}x(λ) := x ∈ X : J(x; λ) ≥ J(x′ ; λ), ∀x′ ∈ X , (8)which corresponds to the set of x attaining the maximum of
the profit function J(·; λ) with a fixed price profile λ. As
shown in [3], x : Rn → Rn in (8) is monotone increasing if
F is convex. Based on this fact, we introduce the following
notion of multidimensional bid functions.
Definition 1: A set-valued function x : Rn → Rn is said
to be a multidimensional bid function with respect to a priceprofile interval [λ] if it is monotone increasing over [λ].

The domain and image of a multidimensional bid function are both multidimensional. This reflects the fact thatprosumption amount transaction is performed at multiple
timeslots. Note that the element with respect to the tthntimeslot, i.e., eTt x : R → R, is not a conventional periodspecific bid function, because it is generally a function of themultidimensional price profile. Note that the graph of eTt x isdepicted as a multidimensional hyperplane. It is not tractablein the conventional bidding system because one-dimensional
bidding curves can only be considered there.

separate. Note, however, that this assumption of separation
does not make sense when an aggregator has energy storage
resources, which provides an ability to shift prosumption
amounts. From this viewpoint, we see that the clearing of
the multiperiod electricity market is not straightforward witha conventional period-specific bidding system. One approachfor systematic market clearing is to devise a method to reducethe multidimensional bid function x in (8) to a separate
approximant. In Section IV-A below, we will develop such an
approximation method for multidimensional bid functions.
III. P ROBLEM F ORMULATION

C. Separability of Multidimensional Bid Functions
For multidimensional bid functions, the following separation will be discussed below.
Definition 2: A multidimensional bid function x : Rn →n
R is said to be separate over a price profile interval [λ] ifeTt x(λ) = xt (λt ),∀t ∈ T(9)

where xt : R → R is a bid function with respect to [λt ].
This notion of separation represents the property that eachelement of a multidimensional bid function is independent
from each other. Thus, a separate multidimensional bid function can be regarded as a simple ensemble of conventional
period-specific bid functions. In addition, another notion of
separation is introduced for convex functions as follows.
Definition 3: A convex function F : Rn → R is said tobe separate over an interval domain [x] if∑
F (x) = t∈T Ft (xt )(10)for an ensemble of Ft : R → R being convex over [xt ].
This corresponds to the additive decomposability of costfunctions. Definitions 2 and 3 have a clear link as follows.
Lemma 1: Consider an aggregator in Section II-A. Then,the multidimensional bid function x : Rn → Rn in (8)is separate over a price profile interval if and only if the
prosumption cost function F in (6) is separate over aninterval domain in the realizable prosumption profile set X
in (4).
Proof: We use the facts from convex analysis theory
shown in Appendix. The definition of the conjugate implies
supx∈X J(x; λ) = F (λ) for J in (7). Furthermore, the first
derivative condition for the supremum leads to λ ∈ ∂F (x).
Because of the equivalence between the conjugates, we seethat x(λ) = ∂F (λ). According to the relation of the subdifferential for separate functions shown in Corollary 31.5.2 of[8], the separation of x over a domain [λ] is rephrased as(∑)x(λ) = ∂λ ∈ [λ](11)t∈T F t (λt ) ,for an ensemble of F t : R → R being convex over [λt ]. Thismeans the separation of F over [λ], which is equivalent tothe separation of F over some [x]. This proves the claim.
Lemma 1 shows that a multidimensional bid function
x is separate over a price interval [λ], i.e., the bidding
curve can be depicted as a one-dimensional curve, only
over a domain [x] where a prosumption cost function F is

A. Remarks on Multiperiod Electricity Market Clearing
In this section, with slight abuse of notation, we denotea symbol of the αth aggregator by that with the subscript
α, e.g., xα and xα . Furthermore, we denote the tuple of asymbol indexed by α ∈ A by that with the subscript A, e.g.,xA := (xα )α∈A . As shown in [3], if at least one prosumption
cost function Fα is smooth, then there exists the unique priceprofile, denoted as λ∗ , such that∏∑∗∃x∗A ∈ α∈A xα (λ∗ ) s.t.
(12)α∈A xα = 0.
In the rest of this paper, we call λ∗ the optimal clearing
price profile and we assume that at least one Fα is smoothfor simplicity.
The optimal clearing of the multiperiod electricity market
can be rephrased as finding λ∗ and x∗A that attain the
prosumption profile balance in (12). This is equivalent tofinding a solution to the convex program of∑∑min α∈A Fα (xα ) s.t.
(13)α∈A xα = 0,
xA

whose Lagrange relaxation is given by{}∑
max min α∈A Fα (xα ) − λT xα .
λ

xA

(14)

As seen here, the Lagrange multiplier appears as the priceprofile λ in (7).
Note that the ISO may be able to solve the convex programof (13) or (14) in a centralized manner. However, this isbased on the premise that every aggregator submits the full
information of the prosumption cost function Fα in (6). Even
though this type of centralized optimization is carried out
in the PJM market [10], such a centralized scheme is notnecessarily desirable from the viewpoint of computational
complexity and privacy of competitive aggregators. Therefore, it is indispensable to design an appropriate bidding
system for the multiperiod electricity market, regarded asa distributed solution scheme for the convex program.
As another approach, dynamic pricing methods can befound in the literature [5], [7]. This approach is mainly
based on the dual ascent algorithm to solve the convex
program in (13), or equivalently (14). Even though such analgorithm can be implemented in a distributed manner, theupdate of each primal variable assumes the strict convexityof the prosumption cost function Fα . In fact, as mentioned
in the end of Section II-A, Fα does not generally satisfy theassumption of the strict convexity.

B. Bidding System Design Problem
For convenience of discussion below, let us introduce thefollowing terminology.
Definition 4: A tuple of prosumption
profiles denoted as∑x̂A is said to be balancing if α∈A x̂α = 0.
In the following, we suppose that each aggregator submits the bidding curves of a separate multidimensional bid
function, i.e., an ensemble of conventional period-specificbidding curves, to the ISO. Under this supposition, theISO can determine a clearing price profile as well as atuple of balancing prosumption profiles. More specifically,
let x̂α : Rn → Rn denote a separate multidimensional bid
function of the αth aggregator, which can be regarded as an
approximant of the original multidimensional bid function
xα given as in (8). Then, there can be uniquely found an
approximate clearing price profile, denoted as λ̂, such that∏x̂A ∈ α∈A x̂α (λ̂)(15)∑
is balancing if α∈A x̂α is strictly monotone increasing.
To evaluate quality of a bidding system, we define thefollowing measure for balancing prosumption profiles.
Definition 5: Let x̂A be a tuple of balancing prosumption
profiles. A deadweight loss with respect to x̂A is defined by∑∆(x̂A ) := α∈A Fα (x̂α ) − F ∗(16)where F ∗ denotes the minimum social cost of the convex
program in (13).
By definition, we clearly see that the deadweight loss isnonnegative and it becomes zero for the optimal solution of
the convex program in (13). On the basis of arguments above,
we formulate a bidding system design problem as follows.
Problem: Consider a set of aggregators in Section II-A.
Devise a market clearing scheme that can find a tuple ofprosumption profiles and a clearing price profile, denoted asx̂A and λ̂, such that the following requirements are satisfied.
– Each aggregator submits an ensemble of bidding curves
associated with a separate multidimensional bid function x̂α to the ISO.
– The ISO determines x̂A and λ̂ such that x̂A is balancing
and (15) holds.
– The deadweight loss with respect to x̂A is small enough.
The major difficulty in this problem is how to determine
a good separate multidimensional bid function x̂α such thatthe resultant deadweight loss is small enough. In the rest of
this paper, we will address this issue.
IV. S EQUENTIAL D ETERMINATION OF S EPARATEM ULTIDIMENSIONAL B ID F UNCTIONS
A. Bid Function Determination at Specific Timeslot
In this subsection, we propose a sequential procedure forderiving a separate multidimensional bid function, denoted
as x̂, which corresponds to an approximant of the original
multidimensional bid function x in (8). For simplicity ofnotation, we drop the subscript α as long as our attention canbe focused on one aggregator without distinction of them.
With the notation of stacked vectors like
xi:j := (xi , xi+1 , . . . , xj )T ,where the subscript is associated with timeslots, we suppose that a price subprofile λ̂1:t−1 as well as a transacted
prosumption subprofile x̂1:t−1 have been determined by the
ISO for the first to (t − 1)th timeslots. Our objective here isto determine a period-specific bid function x̂t : R → R for
the tth timeslot, based on the premise that λ̂1:t−1 and x̂1:t−1
have been determined in advance. Without loss of generality,
they are assumed to be empties when t = 1.
To formulate a profit function with respect to the tth
timeslot, let us consider a cost function in the form of{}ft (xt ) := − min max λT
t+1:n xt+1:n − F (x̂1:t−1 , xt:n )λt+1:n xt+1:n(17)
where the minimization and maximization are subject toλt+1:n ∈ [λ̂t+1:n ],

t−1
∑
τ =1

eτ x̂τ +n∑eτ xτ ∈ X .

(18)

τ =t

In this program, the interval [λ̂t+1:n ], called a price prediction interval, is determined by the aggregator of interest.
Over the feasible domain of ft , denoted by [x̂t ], we consider
a period-specific profit function Jˆt defined as
Jˆt (xt ; λt ) := λt xt − ft (xt ).

(19)Then, we obtain the following period-specific bid function.
Lemma 2: Consider an aggregator in Section II-A anddefine the sequential profit function Jˆt as in (19). For any
price prediction interval [λ̂t+1:n ] and any constants λ̂1:t−1
and x̂1:t−1 , the set-valued function{}x̂t (λt ) := xt : Jˆt (xt ; λt ) ≥ Jˆt (x′t ; λt ), ∀x′t ∈ [x̂t ] (21)is a bid function with respect to any price interval [λt ].
Proof: To prove the claim, we show that x̂t is monotone
increasing. First, we suppose that ft in (17) is a convex
function. Under this supposition, x̂t (λt ) = ∂f t (λt ) holds
where f t denotes the conjugate of ft . Thus, the monotonicityof the subdifferential proves the monotonicity of x̂t .
Let us show the convexity of ft . We notice that{}F t (xt , λt+1:n ) := max λT
t+1:n xt+1:n − F (x̂1:t−1 , xt:n )xt+1:ncorresponds to the partial conjugate of the convex function
F (x̂1:t−1 , xt , ·). As shown in Theorem 33.1 of [8], F t isa saddle function, i.e., it is concave with respect to the firstvariable and convex with respect to the second. Furthermore,
ft is given as the maximum of the collection of convex
functions −F t (·, λt+1:n ), which is shown to be convex; see
Theorem 5.5 of [8]. Thus, ft is convex.
On the basis of the bid function x̂t in (21), each aggregatorcan depict own period-specific bidding curve for the tth
timeslot market. By definition, provided that the feasible
domain [x̂t ] is nonempty, so is [x̂t+1 ] for any x̂t ∈ x̂t (λ̂t )with a clearing price λ̂t ∈ [λt ]. Thus, according to thisB. Multiperiod Electricity Market Clearing
For a set of aggregators, let x̂α : Rn → Rn denote theseparate multidimensional bid function of the αth aggregator
given as in Section IV-A. By aggregating all bidding curves,
the ISO can determine a tuple of prosumption profiles anda clearing price profile, denoted as x̂A and λ̂, such that x̂Ais balancing and (15) holds. This gives one solution to thebidding system design problem in Section III-B.
Let us discuss the resultant deadweight loss, i.e., thequality of the bidding system. In fact, we can guarantee theoptimality at least when every multidimensional bid function
xα is originally separate. This is formally stated as follows.
Theorem 1: Consider a set of aggregators in Section II-A
and let x̂α : Rn → Rn denote the separate multidimensionalbid function of the αth aggregator given as in (21). Assume
that every multidimensional bid function xα in (8) is separate
over a price profile interval [λ]. If λ∗ ∈ [λ] for the optimal
price profile λ∗ in (12), then the clearing price profile λ̂ suchthat x̂A in (15) is balancing is identical to λ∗ . Furthermore,
the deadweight loss with respect to x̂A is zero.
Proof: Let us focus our attention on one aggregator and
we drop the subscript α without distinction of aggregators.
From the assumption that x is separate over [λ], we see that
F is separate over some domain [x] as shown in Lemma 1.
When F is represented as in (10), ft in (17) is identical to
Ft up to an additive constant. Thus, x̂t in (21) is equal to∂F t . This is also equal to the tth element of x represented
as in (11). Therefore, provided that the multidimensional bid
functions of all aggregators are separate and λ∗t ∈ [λt ], theoptimal solutions are to be found, namely λ̂t = λ∗t and x̂t =x∗t , which satisfy x∗t ∈ ∂F t (λ∗t ). This proves the claim.
Theorem 1 can guarantee the optimality of market clearing
under the assumption that every xα is originally separate.
However, such an assumption may not make sense especiallywhen an aggregator has large energy storage resources. This
will be numerically demonstrated in Section IV-C below.
C. Numerical Example
We consider the multiperiod electricity market consistingof three aggregators. On the day of interest, the transacted
prosumption amounts are determined at every 90 minutes,
i.e., the dimension of prosumption profiles is n = 16.
The first and second aggregators are supposed to manage
five thousand residential consumers with energy storage
resources, formulated asxα = −lα + ηαout δαout − η1in δαin ,αα ∈ {1, 2}.

Furthermore, the third∑9aggregator manages nine generators,
formulated as x3 = i=1 gi , where the subscript i represents
the label of generators. In our simulation, we set η1out =η2in = 0.95 and η2out = η2in = 0.94. The inverter and battery
capacities are supposed to be 3.5 [kW] and 10.5 [kWh],× 10 6(a)1.1

Agg. 1
Agg. 2

Cost [JPY]5000
Load [kW]sequential determination of x̂t and λ̂t , each aggregator candetermine a separate approximant x̂ : Rn → Rn as theensemble of x̂t .

4000
3000
2000

0

5

10
15
Time [h]

20

25

(b)

1.05
1
0.95
0.9
0

10

20
30
40
Penetration level [%]50

Fig. 1. (a) Load profiles of aggregators. (b) The resultant social costs
versus levels of energy storage penetration.

respectively, for one residential consumer. The load profiles
l1 and l2 are plotted in Fig. 1(a).
The fuel cost function of thermal generation is supposed
to be a linear function of∑9(22)
G3 (g1 , . . . , g9 ) = i=1 3i × 1T
16 gi ,
which is convex but not strictly convex. The lower and upperbounds of gi are given as
g1 , g2 , g3 ∈ [0, 1500 × 1n ],g7 , g8 , g9 ∈ [0, 800 × 1n ].

g4 , g5 , g6 ∈ [0, 1000 × 1n ],The battery charge and discharge cost functions D1 andD2 are given as follows. Let s0α denote the initial amount of
stored energy, i.e., the state of charge, which is defined as thedeviation from a neutral value. Then, the amount of stored
energy deviation at the termination time is represented asinoutsα (δα ) := s0α + 1Tn (δα − δα ).

(23)

In our simulation, we set s01 = s02 = 0. It is reasonable tosuppose that a higher level of final stored energy is more
preferable than a lower level, and vice versa. To take into
account this aspect, each aggregator is supposed to assess
the value of final stored energy by
Dα (δα ) := −d (sα (δα ))where d : R → R is a concave function given as
a4 (s − s) + a3 s,
s ≤ s,
a3 s,
0 ≤ s < s,d(s) :=as,s≤ s < 0,2
a1 (s − s) + a2 s,
s < s.
We set a1 = 20, a2 = 10, a3 = 6.7, a4 = 3.3, s = −1.31,
and s = 1.31. In this setting, Dα is also convex but notstrictly convex. Then, we obtain a convex prosumption costfunction Fα for each aggregator. We give the price prediction
interval [λ̂2:16 ] in (18) as the multidimensional interval from10 to 23 [JPY/kWh].
Varying the levels of energy storage penetration, we calculate the resultant deadweight loss ∆(x̂A ) in (16) when using
the separate approximants x̂α derived∑as in Section IV-A. In
Fig. 1(b), the resultant social costs α∈A Fα (x̂α ) and F ∗are plotted by the blue solid line with squares and the red
dotted line with circles, respectively. The difference between
them represents the deadweight loss and the horizontal axis
represents the percentage of residential consumers havingenergy storage. From this figure, we see that the deadweightloss is zero at least when there is no energy storage, i.e., the
0% penetration level. This is because x1 , x2 , and x3 in this

15
10
5
0

Fig. 2.

5

10
15Time [h]

20

25

1500
1000
500
0
0

5

10
15
Time [h]

20

25

4000

2000
Stored energy [kWh]Optimal priceApprox. priceCharge/discharge power [kW]Thermal generation [kW]Price [JPY/kWh]20

Agg. 1
Agg. 2

1000
0
-1000
-2000
0

5

10
15
Time [h]

20

25

2000
0
-2000
-4000
0

5

10
15
Time [h]20

25

Price profiles, dispatchable power generation profiles, charge and discharge power profiles, and stored energy profiles in 10% penetration level.

case are all separate, i.e., x̂α is identical to xα as shown in
Theorem 1.
However, we further see that the deadweight loss as well as
the resultant social cost increase as the level of energy storage
penetration increases. This implies that the approximation ofxα by x̂α is not good enough when the level of energy
storage penetration is high. In fact, as shown in Fig. 2, theresultant profiles of clearing prices and decision variables(the blue solid lines) have relatively large discrepancy fromthe optimal ones (the red dotted lines) even in the case of
the 10% penetration level.
V. B IDDING S YSTEM FOR E NERGY S HIFTABILITY
A. Aggregator Model in Multiresolved BasisIn the rest of this paper, for simplicity of discussion, we
suppose that n is a power of 2, i.e., n = 2N for a natural
number N . Under this supposition, every index h ∈ H with
H := {0, 1, . . . , 2N − 1}
can be represented in the binary form of∑(j)
h = j∈J σh 2j , J := {0, 1, . . . , N − 1}(j)where each σh takes a binary value of 0 or 1. As being
compatible with this binary representation, we consider avector of coordinates given as
(0)

(1)

(N −1)

uh := ph ⊗ ph ⊗ · · · ⊗ ph

,

h∈H

(j)
where each ph is the binary vector of

(
)
)
)(
1 (j)
1 (
1
1
(j)
(j)
+ √ σh.
ph := √ 1 − σh
−1
1
2
2

The system of these coordinates forms an orthonormal basis
denoted by {uh }h∈H , which we call a multiresolved basis.
Note that uh with a small number of h corresponds to abasis of low temporal resolution. On the other hand, uh witha large number of h corresponds to a basis of high temporal
resolution. For example,)( )N( )N (
12N −1
1
1
√
√
u0 =
(24)
12N , u1 =2
2
−12N −1
correspond to the zero (lowest) frequency coordinate and thesecond lowest frequency coordinate, respectively. The basis
transformation can be regarded as a bijection between thetime domain and a frequency-like domain.
On the premise of the multiresolved basis, we consider
formulating a transformed multiperiod electricity market
with regard to energy shiftability. To∑ this end, we represent
a prosumption profile as x =h∈H wh uh where whdenotes the hth component in the multiresolved basis. Thisis rewritten as x = U w where w := (wh )h∈H denotes thestacked component vector andU := [u0 , u1 , . . . , u2N −1 ](25)denotes the unitary matrix associated with the multiresolved
basis. With this basis transformation, we can define anNalternative cost function H : R2 → R defined as
H(w) := F (U w).

(26)

This transformed cost function represents the prosumption
cost in the multiresolved basis. For example, its partial
derivative with respect to w0 represents the change rate ofcosts for increasing the total prosumption amount on the
day of interest, and its partial derivative with respect to w1
represents the change rate of costs for shifting prosumption
amounts from the latter half (afternoon) timeslots towards
the former half (morning) timeslots. They correspond to thecoordinates u0 and u1 in (24).
B. Derivation of Separate Multidimensional Bid Functions
Consider the basis transformation of the price profile as
λ = U η, where η denotes the transformed price profile inthe multiresolved basis. Note that each price ηh is associatedwith the corresponding coordinate uh . Then, the transformed
profit function, which is equivalent to J in (7), is defined as
K(w; η) := η T w − H(w).

(27)

On the basis of this profit function, we can define themultidimensional bid function{}w(η) := w : K(w; η) ≥ K(w′ ; η), ∀w′ ∈ U T X , (28)which is given in the same manner as x in (8).
Let us consider deriving a separate approximant of themultidimensional bid function w according to the sequential
procedure in Section IV-A. In the same manner as that for x̂tin (21), we can derive a bid function ŵh : R → R associatedwith the hth resolution market. More specifically, under the
supposition that a price prediction interval [η̂h+1:2N −1 ] isgiven and constants ŵ0:h−1 and η̂0:h−1 are determined, we
define a convex function gh : R → R in the same manner as
ft in (17). Then, for the hth resolution profit functionK̂h (wh ; ηh ) := ηh wh − gh (wh ),(29)

we obtain the bid function as{}ŵh (ηh ) := wh : K̂h (wh ; ηh ) ≥ K̂h (wh′ ; ηh ), ∀wh′ ∈ [ŵh ] ,(30)whose derivation is similar to that of (21).

NN

Finally, we obtain a separate approximant ŵ : R2 → R2
as the ensemble of ŵh . Note that the original multidimensional bid functions x and w are equivalent up to the basis
transformation, while the separate approximants x̂ and ŵ arenot equivalent. Thus, the resultant values of deadweight loss
are different in general.

Furthermore, provided that the price prediction interval[η̂h̄−i+1:h̄ ], which is the domain of the minimization withrespect to ηh̄−i+1:h̄ in (35), involves 0 as a feasible element,
the convex program inside gh̄−i (ŵh̄−i ) has the solution of

C. Energy Shiftability Market Clearing
Consider a set of aggregators. As shown in [3], the optimal
prosumption profile x∗α in (12) of at least one aggregator isshiftable, i.e., the directional derivative [8] of Fα at x∗α withrespect to ei − ej is zero for every pair (i, j) ∈ T × T , if
and only if the optimal price profile λ∗ in (12) lies in theimage of u0 in (24), i.e., λ∗1 = · · · = λn∗ . This implies that the
optimal clearing price profile levels off due to the shiftabilityof prosumption profiles. Note that this price leveling off is
expected when the level of energy storage penetration is
sufficiently high, because energy storage can be seen as aresource having a larger shiftable domain of prosumption
profiles. The following theorem shows that, under sucha situation of price leveling off, the transformed bidding
system in Section V-B can realize the optimal clearing of
the multiperiod electricity market.
Theorem 2: Consider a set of aggregators in Section II-A.
For a price prediction interval such that 0 ∈ [η̂1:2N −1 ] holds,
NN
let ŵα : R2 → R2 denote the separate multidimensionalbid function of the αth aggregator given as in (30). Then,with regard to the clearing price profile η̂ such that∏x̂A = (U ŵα )α∈A , ŵA ∈ α∈A ŵα (η̂)(31)

This is ensured by (34) with θ = i. Therefore, we see thatis balancing, the optimal price profile λ∗ in (12) is given byλ∗ = η̂0 u0 if and only if η̂ satisfiesη̂1 = η̂2 = · · · = η̂2N −1 = 0.

(32)

Furthermore, the deadweight loss with respect to x̂A is zero.
Proof: For simplicity of notation, we drop the subscript
α as long as we can focus on one aggregator. We first provethat (32) is equivalent to
0 ∈ ∂ Ĥh (ŵh ),∀h ∈ {h1 , h2 , . . . , h̄}(33)where h̄ := 2N − 1 and Ĥh (wh ) := H(ŵ0:h−1 , wh , ŵh+1:h̄ ),provided that ŵ ∈ ŵ(η̂) holds as in (31). This proof is doneby induction in the descending order. Initially, we haveK̂h̄ (wh̄ ; ηh̄ ) = wh̄ ηh̄ − H(ŵ0:h̄−1 , wh̄ )for K̂h in (29). Owing to the convexity of H(ŵ0:h̄−1 , ·), thecondition of η̂h̄ = 0 is equivalent to the fact that the h̄thelement ŵh̄ maximizes K̂h̄ (·; 0), namely 0 ∈ ∂ Ĥh (ŵh ) forh = h̄. Next, supposing that0 ∈ ∂ Ĥh (ŵh ),∀h ∈ {h̄ − θ + 1, . . . , h̄}(34)for θ = i, we show that (34) holds again for θ = i + 1. Notethat η̂h̄−i = 0 is equivalent to the fact that ŵh̄−i satisfies
0 ∈ ∂gh̄−i (ŵh̄−i ) where{ T
max ηh̄−i+1:gh̄−i (wh̄−i ) = − minh̄ wh̄−i+1:h̄ηh̄−i+1:h̄ wh̄−i+1:h̄} (35)−H(ŵ0:h̄−i−1 , wh̄−i , wh̄−i+1:h̄ ) .

ηh̄−i+1:h̄ = 0,wh̄−i+1:h̄ = ŵh̄−i+1:h̄ .

0 ∈ ∂gh̄−i (ŵh̄−i ) = ∂ Ĥh (ŵh ),which proves (34) for θ = i + 1. Hence, the equivalence
between (32) and (33) is proven.
Notice that the optimal multidimensional bid function w
in (28) is given as w = ∂H, where H denotes the conjugateof H. For the arguments below, let us confine our attention
on the domain of η ∈ im e0 , where we denote the unit vector
associated with η0 by e0 . Then, for any ŵ0 ∈ eT
0 ∂H(e0 η0 ),it follows that( T)
e0 ∂H(e0 η0 )w(e0 η0 ) =∂H(ŵ0 , 0)
where H(ŵ0 , ·) denotes the conjugate of H(ŵ0 , ·). In thefollowing, we show that()Tw(e0 η0 ) = ŵ0 (η0 ), ŵ1 (0), . . . , ŵh̄ (0) .
(36)Note that the minimization in (35) has the solution being
equal to zero if (32) or equivalently (33) holds. Thus, without
loss of generality, we can assume that [η̂1:h̄ ] = 0 for its
feasible domain. This reduces the formula of K̂h in (29) asK̂h (wh ; ηh ) = wh ηh − max {−H(ŵ0:h−1 , wh:h̄ )}. (37)wh+1:h̄Let us first consider the case of h = 0. From the reduced
formula, we see that ŵ0 (e0 η0 ), which is given as the set of
w0 attaining the maximum of K̂0 (·; e0 η0 ), can be expressedas ŵ0 (η0 ) = eT
0 ∂H(e0 η0 ), which proves the identity of thezeroth element of (36).
On the other hand, because the maximization in (37) has
a solution of ŵh+1:h̄ , ŵh (0) for each h ≥ 1 is given as theset of wh such that 0 ∈ ∂Hh (wh ). This is satisfied for ŵhbecause of (33). Therefore, stacking them for all h ≥ 1, we
have (ŵ1 (0), . . . , ŵh̄ (0))T = ∂H(ŵ0 , 0), which proves theidentity of the remaining elements of (36).
Let us denote the left-hand side of (36) associated with theαth aggregator by wα (e0 η0 ). In this notation, the balance
equation in (31) with (32) implies that
∏
∑
∗
∗
∃wA
∈ α∈A wα (e0 η̂0 ) s.t.
(38)
α∈A wα = 0.
Therefore, from the uniqueness ∑of the optimal price profile,∗
we see that η ∗ satisfying 0 ∈α∈A w α (η ) is given by∗∗η = e0 η̂0 . Hence, λ = η̂0 u0 is proven. The deadweightloss of x̂A is zero because (38) is equivalent to (12).
Theorem 2 implies that, in a situation where the optimal
price profile λ∗ levels off, the optimal price profile isobtained as λ∗1 = · · · = λ∗n = ( √12 )N η̂0 , where η̂0 denotes
the clearing price for the total prosumption on the day. As
discussed in [3], such price levelling off can be expected inhigh penetration of energy storage. From this viewpoint, the

15
10
5
0

Fig. 3.

5

10
15Time [h]

20

25

1500
1000
500
0
0

5

10
15
Time [h]

20

25

4000

2000
Stored energy [kWh]Optimal priceApprox. priceCharge/discharge power [kW]Thermal generation [kW]Price [JPY/kWh]20

Agg. 1
Agg. 2

1000
0
-1000
-2000
0

5

10
15
Time [h]

20

25

2000
0
-2000
-4000
0

5

10
15
Time [h]20

25

Price profiles, dispatchable power generation profiles, charge and discharge power profiles, and stored energy profiles in 10% penetration level.

sequential determination of separate multiperiod bid functions in the multiresolved basis can work well for clearing
the multiperiod electricity market under high energy storage
penetration. Note that the sequential determination does notwork well in the original time basis as demonstrated in
Section IV-C.
D. Numerical Example
Consider the same setting as that in Section IV-C. We givethe price prediction interval [η̂1:15 ] as the multidimensional
interval from −9.5 to 9.1. Varying the levels of energy
storage penetration, we calculate the resultant deadweightloss for the transformed ∑bidding system in Section V-B.
The resultant social cost α∈A Hα (ŵα ) is over plotted inFig. 1(b) by the green solid line with diamonds. From thisfigure, we see that the resultant social cost decreases as
the penetration level increases. Furthermore, the deadweight
loss, represented as the difference from the red dotted linewith circles, is smaller than the deadweight loss in Section IV-C. In fact, in the case of no energy storage penetration, i.e., the 0% penetration level, the resultant deadweightloss is shown to be zero because the separate approximantŵα is identical to the original multidimensional bid function
wα for every α ∈ {1, 2, 3}. In addition, as expected from
Theorem 2, the deadweight loss becomes smaller as thepenetration level becomes higher. For reference, the resultant
profiles of clearing prices and decision variables in the caseof the 10% penetration level are plotted in Fig. 3, where weuse the same legends as those in Fig 2. From these results,
we see that the distributed approximate scheme for market
clearing has good compatibility with the bidding system inthe multiresolved basis.
VI. C ONCLUDING R EMARKSIn this paper, we have designed a bidding system for a
multiperiod electricity market with consideration of the pricing of the shiftability of energy storage. More specifically,
we have first developed a distributed approximate scheme
for multiperiod electricity market clearing. Then, based ona basis transformation similar to the Fourier transformation,
we propose a bidding system with explicit consideration of
the pricing of energy shiftability. It has been theoreticallyand numerically shown that the distributed approximate
scheme has good compatibility with the bidding systemin the Fourier-like basis. Generalization to the integrationof renewable power generation, which may have a steep
fluctuation, is a meaningful future work to pursue.

R EFERENCES[1] A. M. Annaswamy and M. Amin, “IEEE vision for smart grid controls:
2030 and beyond,” 2013.
[2] B. A. Hamilton, J. Miller, and B. Renz, “Understanding the benefitsof the smart grid-smart grid implementation strategy,” United States:
United States Department of Energy s National Energy Technology
Laboratory, 2010.
[3] T. Ishizaki, M. Koike, N. Yamaguchi, A. Ueda, T. Bo, and J.-i. Imura,“A distributed scheme for power profile market clearing under highbattery penetration,” in Proceedings of the 19th IFAC World Congress(to appear), 2017.
[4] G. He, Q. Chen, C. Kang, P. Pinson, and Q. Xia, “Optimal bidding
strategy of battery storage in power markets considering performancebased regulation and battery cycle life,” 2015.
[5] J. Hansen, J. Knudsen, and A. Annaswamy, “A dynamic market
mechanism for integration of renewables and demand response,” IEEE
Tran. on Smart Grid, 2015.
[6] G. Liu, Y. Xu, and K. Tomsovic, “Bidding strategy for microgrid
in day-ahead market based on hybrid stochastic/robust optimization,”
Smart Grid, IEEE Transactions on, vol. 7, no. 1, pp. 227–237, 2016.
[7] D. J. Shiltz, M. Cvetković, and A. M. Annaswamy, “An integrated
dynamic market mechanism for real-time markets and frequency
regulation,” IEEE Transactions on Sustainable Energy, vol. 7, no. 2,
pp. 875–885, 2016.
[8] R. T. Rockafellar, Convex analysis. Princeton university press, 1970.
[9] D. T. Luc, “On generalized convex nonsmooth functions,” Bulletinof the Australian Mathematical Society, vol. 49, no. 1, pp. 139–149,
1994.
[10] A. L. Ott, “Experience with PJM market operation, system design,
and implementation,” IEEE Transactions on Power Systems, vol. 18,
no. 2, pp. 528–534, 2003.

A PPENDIX
We overview several facts from convex analysis theory [8].
The conjugate of a convex function F is defined by{ T}
λ x − F (x) ,(39)
F (λ) := supx∈dom F

where dom F denotes the effective domain of F . It is known
that F is convex and the conjugate of F coincides with F
as long as F is convex. Furthermore, the strict convexity of
F is equivalent to the smoothness of F . The transformationin (39) is called the Legendre-Fenchel transformation.
The subdifferential of a convex function F is defined by∂F (x) := {λ : F (x′ ) ≥ F (x) + λT (x′ − x), ∀x′ ∈ dom F },which is a set-valued function with a convex image. Corollary 31.5.2 of [8] shows that ∂F : Rn → Rn is monotone
increasing if F is convex. In particular, it is shown to bestrictly monotone increasing if F is strictly convex; see
Theorem 2.1 in [9]. Furthermore, x ∈ ∂F (λ) and λ ∈ ∂F (x)are equivalent as shown in Theorem 23.5 in [8].