fix step size 過(guò)高引起的Precision異常
Projection method 里面分為兩步嘁信,一步是步長(zhǎng)測(cè)試尤蛮,二是其他函數(shù)值得計(jì)算呐馆。從fix step size做起狞尔,可以先排除step size adjust 上的問(wèn)題丛版。
數(shù)據(jù)準(zhǔn)備好,初次實(shí)驗(yàn)成功偏序,下面看邊界條件页畦。
回想Fix 遇到什么問(wèn)題?步長(zhǎng)超過(guò)邊界不能迭代
Ptt并不高
step size 沒(méi)有下限研儒,怎么低都可以豫缨,就是計(jì)算效率不高。比如極限情況殉摔,步長(zhǎng)為0.0001州胳,迭代數(shù)百萬(wàn)次记焊,精度一直在提高逸月。迭代數(shù)百萬(wàn)次,說(shuō)明運(yùn)行穩(wěn)定遍膜。精度一直在提高說(shuō)明碗硬,雖然慢瓤湘,但沒(méi)有錯(cuò)誤。
下面看看另一個(gè)極端情況恩尾。
發(fā)現(xiàn)Bug:精度突變弛说。最后出現(xiàn)古怪符號(hào)。精度最高e-16
繼續(xù)測(cè)試步長(zhǎng)極限
- Note: step size=40, Exception appears, the precision reduce and then, exception appears.
- Note: step size=50, Exception again appears, the precision reduce and then, exception appears.
-
note: we have not good results when step size or precision go to limitations.
Route cost 沒(méi)有問(wèn)題
Mode cost 在前面的測(cè)試中都是正常的翰意, 這里也沒(méi)有出現(xiàn)異常木人。只是數(shù)值比path cost要大一些。
從程序上看不出問(wèn)題
測(cè)試flow update冀偶,從結(jié)果上看不出問(wèn)題醒第,就是一個(gè)projection 操作,從中間過(guò)程看进鸠,也沒(méi)有問(wèn)題稠曼。首先是mode做流量調(diào)整,fix step size 乘以 transfer cost客年。投影霞幅。流量守恒。然后解路徑流量量瓜。將升級(jí)后的流量分配到下轄的各條路徑上司恳。按mode demand plus 守恒。投影運(yùn)算榔至,流量-cost抵赢,沒(méi)錯(cuò)。
然后計(jì)算誤差
誤差沒(méi)有縮小唧取,反而變大了铅鲤。但是結(jié)合以前的情況也沒(méi)有問(wèn)題。在接下來(lái)的幾次迭代中枫弟,精度反復(fù)波動(dòng)邢享。第七次迭代,誤差急劇變大淡诗。
第八次:not a number
path flow -- path cost (OK) -- mode cost (Exception) ---
檢查mode cost
nest sum 出現(xiàn)有的nest一點(diǎn)流量都沒(méi)分到的情況骇塘。
在上一步中所有流量都分到一個(gè)nest mode上。上一步就有問(wèn)題了韩容,流量分配不對(duì)款违。logit based 每種方式都有一定的流量。
ptt (OK)
mode cost :
input是mode demand群凶,其中一個(gè)demand是0插爹,導(dǎo)致endogenous cost為負(fù)無(wú)窮。取值-999999999999999。外部cost是0.603289. 所以赠尾,總體cost是:
然后力穗,減去最小值就變成
這個(gè)基本是沒(méi)錯(cuò)的。有點(diǎn)誤差气嫁。不知道為啥当窗,小數(shù)點(diǎn)后面的問(wèn)題。簡(jiǎn)單的減法寸宵,小數(shù)點(diǎn)后面幾位不一致崖面。
正常情況下沒(méi)有這個(gè)問(wèn)題,可能是數(shù)值太大引起的梯影∷恢欤可以單獨(dú)領(lǐng)出來(lái)試試。
總的來(lái)說(shuō)就是Demand=0引起的光酣。
那么就是怎么處理demand=0的情況?
在輸入端控制救军,給一個(gè)極小的demand:e-20财异。
if (nm[nestmode_i].NestSum != 0) {
nm[nestmode_i].EndogenousCost = \
(MuE / Theta)*log(0.00000000000000000001/ pow(nm[nestmode_i].Membership, 1 / MuE)) \
+ ((1 - MuE) / Theta)*log(nm[nestmode_i].NestSum);
}
else {
nm[nestmode_i].EndogenousCost = \
(MuE / Theta)*log(0.00000000000000000001 / pow(nm[nestmode_i].Membership, 1 / MuE)) \
+ ((1 - MuE) / Theta)*log(0.00000000000000000002);
}
Demand 調(diào)整以后,再測(cè)試精度唱遭。在第八次迭代出現(xiàn)not a number的情況戳寸。
在求ptt和path cost的時(shí)候應(yīng)該沒(méi)有問(wèn)題了。
問(wèn)題可能處在termination test上拷泽。
IndexNM[nestmode_i] = IndexNM[nestmode_i] + (p[path_i].Pf / nm[nestmode_i].Demand) * (p[path_i].PttTrans / p[path_i].Ptt);
分母nm[nestmode_i].Demand)確實(shí)可能為0疫鹊;
至此,基本確定了fix step size 過(guò)高引起的Precision異常的問(wèn)題司致。
Self-adaptive Projection 的問(wèn)題拆吆,在ND上還是有問(wèn)題,(Demand)
FW 解ND net 已經(jīng)可以了脂矫。
Fix step size 至少處理了
S2FunValueNth(l, p, nm, od);
S3SolutionUpdate(p, nm, od, Alpha);
TerminationIndex = S8Termination(p, nm, od); printf("%.10lf\n", TerminationIndex);
S9Transit(p, nm);
都沒(méi)什么問(wèn)題枣耀,其中S2和S3是重要步驟,S8和S9相對(duì)簡(jiǎn)單庭再。
在self adaptive中捞奕,還有更多的計(jì)算量。
2拄轻,3颅围,4,5恨搓,6試算步長(zhǎng)
7調(diào)整步長(zhǎng)之Gamma調(diào)整
8收斂測(cè)試
9流量交替
9.1步長(zhǎng)交替
先考慮那些東西要交替院促?flow/demand;alpha ;Gamma每次賦值就行一疯。flow/demand cost都是計(jì)算而來(lái)的。在SAGP中涉及到的所有變量夺姑,x, y, alpha, alpha Plus, 這些變量中需要進(jìn)行交替操作的只有input(流量)和步長(zhǎng)墩邀。這兩步驟很簡(jiǎn)單,且經(jīng)過(guò)多次操作沒(méi)有問(wèn)題盏浙。
8收斂測(cè)試剛剛做過(guò)眉睹,沒(méi)問(wèn)題。
7步長(zhǎng)調(diào)整中Gamma步長(zhǎng)是步長(zhǎng)調(diào)整階段的變量废膘。其他還有一個(gè)變量是最小整數(shù)l竹海。S2-S6就是為了找到最優(yōu)l,S7就是為了調(diào)整這個(gè)Gamma
S7和S6幾乎相同丐黄,只是S6中的Delta是外部決定的斋配,在S7直接給了0.5. 2-Delta始終大于0.5,所以S6中的不等式比S7更容易達(dá)到灌闺。如果S7中的不等式也成立艰争,說(shuō)明這個(gè)不等式條件很容易達(dá)成,步長(zhǎng)可以擴(kuò)大一點(diǎn)桂对。由于S7和S6幾乎相同甩卓,S7的檢查也pass。
檢查S6. 首先大于等于蕉斜,基本公式都對(duì)過(guò)逾柿,沒(méi)有問(wèn)題。
S6中有三大函數(shù):都是內(nèi)積宅此。內(nèi)積的input分別是demand/flow机错;demand/flow cost;demand/flow cost (alpha)父腕。先就這內(nèi)積的形式來(lái)看看毡熏。三個(gè)內(nèi)積的計(jì)算形式都很簡(jiǎn)單。沒(méi)有問(wèn)題侣诵,下面就是看看input對(duì)不對(duì)痢法。
S5算demand/flow cost plus(alpha)
S4算demand/flow cost plus
S3算demand/flow plus
S2算demand/flow cost
S3 and S2 已經(jīng)完成測(cè)試,沒(méi)有太多問(wèn)題杜顺。
S4 S5在一些情況下也正常工作财搁,下面就極端情況做個(gè)別測(cè)試即可。FW就是在這種情況下測(cè)出問(wèn)題的躬络。
嘗試跑起來(lái)
initial step size
when initial step size =5尖奔; step size does not change, but after 395 iterations, the precision collapse and go into endless loop
when initial step size =10, the step size keeps increase, but after 94 iterations, the precision collapse and step size becomes zero, but didn't go into endless loop
check above phenomenon
in 93th iteration, we obtain good results, and then these results will be like input.
the input is x, alpha k, gamma, the following y, alpha k plus,x plus, residual is the product of x, alpha k, and gamma.
so, the red part is the input variable, plus gamma and alpha k
-
we can see, Gamma=alpha k+1 / 0.9. match the definition.
as definition, Part1 - Part2 >= Part3 holds.
to be more detail, it can be observed that part 1,2 and 3 are all small.
- input is mode demand/path flow, the output is path flow cost.
path flow ---link flow ---link cost---path cost---- under each nest mode (path cost - min path cost) - input is min path cost under each mode, the output is nest mode cost
min path cost under each mode + nest mode demand + dispersion parameter: Mu, Theta, Membership---nest mode cost
the input is OK, nest mode demand is not equal to 0; so, the input is ok.
note: test whether we can use scientific notation for calculation. Yes, scientific notation can be used for calculation. --> mode cost, mode cost plus, we all use scientific notation to indicate a small enough value to replace zero. - solution update:
-
input: x, y, alpha k+1
demand update is normal, flow conservation is held.
- flow demand
- y plus
input is x plus,
path flow plus---link flow ---link cost ---path cost plus---- under each nest mode (path cost plus - min path cost plus)
min path cost under each mode plus + nest mode demand plus + dispersion parameter: Mu, Theta, Membership---nest mode cost plus
目前為止都沒(méi)問(wèn)題。總結(jié):Debug要回溯提茁,逆序淹禾,從問(wèn)題查起。 -
start from the problem.
only one modification was made, give demand which supposed to be zero to be a very small number.
in 95th iteration, the precision do not collapse, but the step size become zero.
where is the problem? Input or operation?
-
see the operation in S6
part 1 =0;
part 2 ~~0;
part 3 =0;
before 95th iteration, part 1 !=0;
so, check part 1.
It is observed that cost p is unnormal. As expected, cost P in the previous iteration is stable.
S4FunValueNPth(l, p, nm, od); produce unnormal cost p
為什么會(huì)有兩個(gè)相同的最小值
Demand沒(méi)有變異铃岔,正常,說(shuō)明輸入正常峭火。
第一次接待就有問(wèn)題毁习,所以看第一次就好。
如果第一個(gè)mode 的cost加上一個(gè)極小值衅枫,效果如何恐仑?
nm[0].Cost = nm[0].Cost + 1e-100;
1e-100太小,沒(méi)有效果为鳄。這步就先不處理歧斟,其實(shí)本質(zhì)上是精度處理能力不足引起的。
那解決的辦法一般是用科學(xué)計(jì)數(shù)法來(lái)算了偏形;
對(duì)于精度問(wèn)題静袖,想起前幾天一件事,18個(gè)9減一個(gè)極小值俊扭,如0.15732队橙,沒(méi)有得出正確結(jié)果。
為什么cost P變異
bug出現(xiàn):存在多條最短路的時(shí)候花枫,把非負(fù)最短路的流量加起來(lái)刻盐,就不對(duì)。只能選一條最短路劳翰。如果兩條最短路都考慮敦锌,那么做流量升級(jí)的時(shí)候憋肖,沒(méi)有好的升級(jí)方法來(lái)使得流量微調(diào)玛迄。所以只能選一條最短路。
解決辦法呵萨,增加一個(gè)指標(biāo)作為最短路指示指標(biāo)溺蕉。要做任何改動(dòng)必須改一步試一步。
其實(shí)同理,在path flow update上也一樣朽色。
精度是一個(gè)反復(fù)升降的過(guò)程
精度的升降與是否存在多條最短路沒(méi)有關(guān)系
總之,精度是一個(gè)反復(fù)升降的過(guò)程可能是在精度達(dá)到瓶頸以后遇到的問(wèn)題梢褐,這里先不管旺遮。
精度最高到e-20
appendix
debug s2 s3
for (path_i = 0; path_i < NoPs; path_i++) {
//printf("%d %d %d %d %lf %lf\t%lf\t%lf\n", p[path_i].Path[0], p[path_i].Path[5], p[path_i].Nest, p[path_i].Mode, \
p[path_i].Pf, p[path_i].Ptt, p[path_i].PttTrans);
}
S22FunValueNth_ModeFun(nm, od);
for (nestmode_i = 0; nestmode_i < NoNMs; nestmode_i++) {
//printf("%d %d %d %d %lf %lf %lf %lf\n", nm[nestmode_i].O, nm[nestmode_i].D, \
nm[nestmode_i].Nest, nm[nestmode_i].Mode, \
nm[nestmode_i].Demand, nm[nestmode_i].Cost, nm[nestmode_i].EndogenousCost, nm[nestmode_i].ExogenouCost);
}
//S3SolutionUpdate(p, nm, od, Alpha);
for (od_i = 0; od_i < NoODs; od_i++) {
od[od_i].NMNonMinSum2 = 0;
od[od_i].NMNonMinSum = 0;
od[od_i].NMCounter = 0;
od[od_i].Flag = 0;
}
for (nestmode_i = 0; nestmode_i < NoNMs; nestmode_i++) {
if (nm[nestmode_i].CostTrans != 0) {
nm[nestmode_i].DemandP = nm[nestmode_i].Demand - Alpha[1] * nm[nestmode_i].CostTrans;
if (nm[nestmode_i].DemandP < 0) {
nm[nestmode_i].DemandP = 0;
}
//printf("%d %d %d %d %.5e\t%.5e\t%.5e\n", nm[nestmode_i].O, nm[nestmode_i].D, nm[nestmode_i].Nest, nm[nestmode_i].Mode\
, nm[nestmode_i].DemandP, nm[nestmode_i].Demand, Alpha[1] * nm[nestmode_i].CostTrans);
}
}
for (od_i = 0; od_i < NoODs; od_i++) {
for (nestmode_i = 0; nestmode_i < NoNMs; nestmode_i++) {
if (nm[nestmode_i].CostTrans != 0 \
&& od[od_i].O == nm[nestmode_i].O && od[od_i].D == nm[nestmode_i].D) {
od[od_i].NMNonMinSum = od[od_i].NMNonMinSum + nm[nestmode_i].DemandP;
}
}
}
for (od_i = 0; od_i < NoODs; od_i++) {
for (nestmode_i = 0; nestmode_i < NoNMs; nestmode_i++) {
if (nm[nestmode_i].CostTrans == 0 && od[od_i].NMCounter == 0 && \
od[od_i].O == nm[nestmode_i].O && od[od_i].D == nm[nestmode_i].D) {
nm[nestmode_i].DemandP = (od[od_i].Dem - od[od_i].NMNonMinSum);
od[od_i].NMCounter = 1;
}
}
}
for (nestmode_i = 0; nestmode_i < NoNMs; nestmode_i++) {
//printf("%d %d %d %d %.5e\n", nm[nestmode_i].O, nm[nestmode_i].D, nm[nestmode_i].Nest, nm[nestmode_i].Mode, nm[nestmode_i].DemandP);
}