求解Fisher市場(chǎng)均衡問(wèn)題的內(nèi)點(diǎn)算法

畢紅梅， 劉妙華， 趙學(xué)軍

(空軍工程大學(xué)基礎(chǔ)部，西安，710051)

早些時(shí)候，Ye提出了求解Fisher市場(chǎng)均衡問(wèn)題的原-對(duì)偶內(nèi)點(diǎn)算法[1]。但該算法跟蹤的加權(quán)中心路徑并不連續(xù)，從初始可行點(diǎn)出發(fā)還需要通過(guò)勢(shì)函數(shù)下降算法才能尋找到滿足鄰域要求的可行點(diǎn)。近年來(lái)，Potra將Fisher市場(chǎng)均衡問(wèn)題抽象為更一般的線性權(quán)互補(bǔ)問(wèn)題，根據(jù)已知可行點(diǎn)構(gòu)造了一條光滑的中心路徑，給出了求解線性權(quán)互補(bǔ)問(wèn)題的內(nèi)點(diǎn)算法[2]。當(dāng)權(quán)向量為0時(shí)，權(quán)互補(bǔ)問(wèn)題退化為一般的互補(bǔ)問(wèn)題；對(duì)于權(quán)向量部分為0、部分大于0的情形，算法設(shè)計(jì)和分析比一般的互補(bǔ)問(wèn)題更為復(fù)雜。自此以后，求解權(quán)互補(bǔ)問(wèn)題逐漸成為優(yōu)化領(lǐng)域研究熱點(diǎn)之一。文獻(xiàn)[3～4]設(shè)計(jì)了光滑牛頓算法求解權(quán)互補(bǔ)問(wèn)題。文獻(xiàn)[5]在一定假設(shè)條件下提出了線性權(quán)互補(bǔ)問(wèn)題寬鄰域的路徑跟蹤算法。文獻(xiàn)[6～8]提出了求解線性權(quán)互補(bǔ)問(wèn)題的全牛頓步內(nèi)點(diǎn)算法。線性權(quán)互補(bǔ)問(wèn)題來(lái)源于Fisher市場(chǎng)均衡模型，國(guó)內(nèi)外文獻(xiàn)重點(diǎn)討論了線性權(quán)互補(bǔ)問(wèn)題的求解方法和理論證明，算法測(cè)試大多自行設(shè)計(jì)，有些算法的設(shè)定條件并不適用于解決Fisher市場(chǎng)均衡這一實(shí)際問(wèn)題，一些初步的研究結(jié)果可參考文獻(xiàn)[9]。

1 Fisher市場(chǎng)均衡問(wèn)題

在 Fisher 模型中，市場(chǎng)中包含消費(fèi)者C和生產(chǎn)者P。消費(fèi)者i∈C有預(yù)算wi>0從生產(chǎn)者手中購(gòu)買商品使得個(gè)人效用函數(shù)最大化。價(jià)格均衡即給商品定價(jià)，使得每個(gè)消費(fèi)者買到最多的商品，同時(shí)市場(chǎng)出清，即所有的預(yù)算消費(fèi)完并且所有的商品銷售完。不失一般性，假定生產(chǎn)者j有1單位的某種商品出售。

Eisenberg和 Gale證明了市場(chǎng)出清價(jià)格為下述優(yōu)化問(wèn)題的最優(yōu)拉格朗日乘子[10]：

(1)

式中：uij≥0為給定的消費(fèi)者i對(duì)生產(chǎn)者j生產(chǎn)商品的效應(yīng)系數(shù)；xij為消費(fèi)者i從生產(chǎn)者j手中購(gòu)買的商品數(shù)量。

更一般地，式(1)可以改寫成：

(2)

x≥0

假定最優(yōu)解集非空，則式(2)的KKT 條件可以寫成線性權(quán)互補(bǔ)問(wèn)題：

Ax=b,x≥0

-ATy+s=0,s≥0

(3)

xs=w

式中:y、s分別為拉格朗日乘子。

嚴(yán)格可行域定義為：

F++={(x,s)>0:Ax=b,s=ATy}

給定初始嚴(yán)格可行點(diǎn)(x0,y0,s0)∈F++，令c=x0s0，構(gòu)造

w(t)=(1-t)w+tc,t∈(0,1]

(4)

文獻(xiàn)[2]將式(4)替換(3)中的第3個(gè)式子得到新的中心路徑：

Ax=b,x≥0
-ATy+s=0,s≥0
xs=w(t)

(5)

當(dāng)t→0時(shí)，w(t)→w得到(3)的解。

2 原-對(duì)偶內(nèi)點(diǎn)算法

2.1 搜索方向

令γ=min (c),α=βγ, 0<β<1，定義鄰域:N(α,t)={(x,y,s)∈F++:‖xs-w(t)‖≤αt}。

令t+=(1-θ)t，θ∈(0,1)。搜索方向由兩個(gè)方向構(gòu)成，其中仿射方向(Δax,Δay,Δas)由下面的系統(tǒng)給出：

AΔax=0

-ATΔay+Δas=0

(6)

sΔax+xΔas=w-xs

中心方向(Δcx,Δcy,Δcs)由下述系統(tǒng)給出：

AΔcx=0

-ATΔcy+Δcs=0

(7)

sΔcx+xΔcs=c-xs

仿射方向作為預(yù)測(cè)方向，中心方向作為矯正方向。本文采用的中心方向把迭代點(diǎn)拉向c=x0s0以保持可行性，而文獻(xiàn)[2]中的中心方向希望迭代點(diǎn)向w(t)靠近。通過(guò)二分法搜索滿足條件的最大θ值使得新的迭代點(diǎn):

(x+,y+,s+)=(x,y,s)+(1-t+)(Δax,Δay,Δas)+

t+(Δcx,Δcy,Δcs)

(8)

屬于新的鄰域:N(α,t+)={(x+,y+,s+)∈F++:‖x+s+-w(t+)‖≤αt+}。

2.2 算法

算法1 線性權(quán)互補(bǔ)問(wèn)題內(nèi)點(diǎn)算法。

1)輸入ε>0，θ∈(0，1)，t0=1，初始可行點(diǎn)(x0,y0,s0)∈N(α,t0)，令k=0。

2)若‖xksk-w‖<ε，算法停止；否則求解式(6)、(7)得搜索方向:(Δaxk,Δayk,Δask)及(Δcxk,Δcyk,Δcsk)。

3)搜索滿足鄰域條件的θ值θk，令:tk+1=(1-θk)tk;xk+1=xk+(1-tk+1)Δaxk+tk+1Δcxk；yk+1=yk+(1-tk+1)Δayk+tk+1Δcyk；sk+1=sk+(1-tk+1)Δask+tk+1Δcsk。

4)令k=k+1，轉(zhuǎn)2)。

2.3 算法收斂性分析

算法可行需要保證新的迭代點(diǎn)在新的鄰域內(nèi)，并且隨著t的減少，x+s+與w的距離不斷下降，最終滿足精度要求。

假定當(dāng)前迭代點(diǎn)嚴(yán)格可行且滿足‖xs-w(t)‖≤αt。

引理1xs≥(γ-α)te,其中e為單位向量。

證明 由xs≥w(t)-αte≥tc-αte≥(γ-α)te，即得。令：

Δx=(1-t+)Δax+t+Δcx,Δs=(1-t+)Δas+t+Δcs。則x+=x+Δx,s+=s+Δs,于是:

sΔx+xΔs=(1-t+)(w-xs)+t+(c-xs)=

(1-t+)w+t+c-xs=w(t+)-xs

(9)

進(jìn)而：

x+s+-w(t+)=(x+Δx)(s+Δs)-w(t+)=

xs+sΔx+xΔs+ΔxΔs-w(t+)=ΔxΔs

(10)

引理2[11]設(shè)u,v∈Rn,uTv≥0。 記U=diag(u),V=diag(v)，則‖UVe‖≤2-3/2‖u+v‖2。

證明式(9)兩邊同乘以(xs)-1/2，得:

x-1/2s1/2Δx+x1/2s-1/2Δs=(xs)-1/2(w(t+)-xs)。

又：‖w(t+)-xs‖≤‖w(t+)-w(t)+w(t)-xs‖≤

‖w(t+)-w(t)‖+‖w(t)-xs‖≤θt‖(w-c)‖+tα。

再由引理1即得。

要使‖x+s+-w(t+)‖≤αt+=(1-θ)tα，只需

(11)

α=βγ,令β=2/3，則：

假定θ‖w-c‖≤α/8，由式(11)可計(jì)算出θ的下界為

定理1給定ε>0,初始點(diǎn)(x0,y0,s0)∈N(α,t0)。若系統(tǒng)(3)存在最優(yōu)解，則算法1至多經(jīng)過(guò)

證明x+s+與w的距離

‖x+s+-w‖=‖x+s+-w(t+)+w(t+)-w‖≤

‖x+s+-w(t+)‖+‖w(t+)-w‖≤(1-θ)αt+(1-θ)t‖w-c‖=(1-θ)t(α+‖w-c‖)≤(1-θ0)t(α+‖w-c‖)=(1-θ0)k(α+‖w-c‖)。

要使‖x+s+-w‖≤ε，只需(1-θ0)k(α+‖w-c‖)≤ε。

兩邊取對(duì)數(shù)得到:klog (1-θ0)+log (α+‖w-c‖)≤logε。

3 數(shù)值實(shí)驗(yàn)

本節(jié)用算法1求解Fisher 市場(chǎng)均衡問(wèn)題。假定市場(chǎng)上有nc≥2個(gè)消費(fèi)者，np≥2個(gè)生產(chǎn)者，為簡(jiǎn)單起見(jiàn)，不妨設(shè)nc=np=p。記e=ones(1,p)，

對(duì)比式(1)、(2)可知

x=(x11,…,x1p,…,xp1,…,xpp,u1,…,up)T，

b=(e,01×p)T，其中01×p為1行p列0向量，

A為m×n=2p×(p2+p)矩陣

w=(01×p2,w1,…,wp)，wi,uij,i,j=1,2,…，p在區(qū)間(0，1)中隨機(jī)生成，初始點(diǎn)(x0,y0,s0)的取法與參考文獻(xiàn)[1]一致。

用MATLAB R2010b編程實(shí)現(xiàn)算法, 用Intel Xeon2 CPU(3.3 GHz), 8 GiB ram的微機(jī)測(cè)試。當(dāng)滿足條件‖xs-w‖≤10-5時(shí),算法終止。

為說(shuō)明Fisher市場(chǎng)均衡問(wèn)題求解過(guò)程，以最簡(jiǎn)單的p=2時(shí)情形為例，即市場(chǎng)上僅有2位消費(fèi)者和2位生產(chǎn)者。系數(shù)矩陣A為m×n=4×6矩陣

b=[1,1,0,0],權(quán)向量即預(yù)算w=[0,0,0,0,0.957 2,0.485 4],

即2位消費(fèi)者預(yù)算分別為0.957 2和0.485 4，初始點(diǎn)為x0=[0.5,0.5,0.5,0.5,0.471 1,0.668 7]，y0=[2.871 5,2.871 5,1.523 9,1.073 5]，s0=[1.652 0,2.655 3,2.418 8,1.888 5,1.523 9,1.073 5]。

用算法1求解上述問(wèn)題，迭代8次，耗時(shí)0.001 3 s后得到最優(yōu)解x*=[1,0,0,1,0.800 3,0.915 7]，y*=[0.957 2,0.485 3,1.196 0,0.530 0]，s*=[0,0.787 5,0.261 8,0,1.196 0,0.530 0]。

結(jié)果表明消費(fèi)者1購(gòu)買了生產(chǎn)者1生產(chǎn)的1個(gè)單位的商品，消費(fèi)者2購(gòu)買了生產(chǎn)者2生產(chǎn)的1個(gè)單位的商品，相應(yīng)的效應(yīng)系數(shù)分別為0.800 3，0.915 7，若商品價(jià)格定為1.196，0.53，對(duì)應(yīng)相乘得到購(gòu)買商品消費(fèi)金額分別為0.957 2，0.485 4，與預(yù)算一致，商品售完，市場(chǎng)出清。

下面對(duì)每種規(guī)模產(chǎn)生10個(gè)問(wèn)題，計(jì)算10次結(jié)果的平均迭代次數(shù)和平均時(shí)間。表1和表2分別給出了算法1和Ye[1]、Potra[2]算法的迭代步數(shù)和迭代時(shí)間。

表1 Fisher均衡問(wèn)題Ye、Potra及算法1迭代步數(shù)比較

表2 Fisher均衡問(wèn)題Ye、Potra及算法1迭代時(shí)間比較

由表1和表2可見(jiàn)，算法1與 Potra算法在迭代次數(shù)及時(shí)間上都優(yōu)于Ye 算法；當(dāng)問(wèn)題規(guī)模達(dá)到m×n=50×650時(shí)，算法1略優(yōu)于 Potra算法。

4 結(jié)語(yǔ)

本文從求解Fisher市場(chǎng)均衡這一實(shí)際問(wèn)題出發(fā)，設(shè)計(jì)線性權(quán)互補(bǔ)問(wèn)題的內(nèi)點(diǎn)算法，就中小規(guī)模問(wèn)題給出了數(shù)值結(jié)果。對(duì)于大規(guī)模問(wèn)題，迭代步數(shù)和計(jì)算時(shí)間會(huì)有所增加，因此如何調(diào)節(jié)仿射方向和中心方向參數(shù)以提高算法效率，或設(shè)計(jì)更為高效的寬鄰域內(nèi)點(diǎn)算法是未來(lái)的研究重點(diǎn)。