最近根据张欣老师编著的《可计算一般均衡模型的基本原理与编程》学习CGE模型。中间遇到一个疑问，希望大家能探讨下。
书中第P106-107页，是有CES嵌套函数的CGE模型编程（程序语言附后）。在该模型中，资本（QKS）、劳动力（QLS）作为外生变量。
在模拟外界冲击的时候，劳动力增长3倍、资本增长到9倍（QLS.fx=3*QLS0; QKS.fx=9*QKS0;）就不能得到均衡结果。出现错误‘7 error(s): vcPower: FUNC DOMAIN: x**c, x < 0 (RETURNED 0)’。这是什么原因呢，该如何解决呢？

我自己根据中国2007年的投入产出表做的静态CGE模型，资本冲击达到3倍以上就无法得到均衡解。改成动态CGE能解决这个问题吗？
希望有经验的老师、同学帮忙解答。非常感谢！

张老师书中的程序语言：

set ac /agri, manu, serv, lab, cap, hh, total/;
set a(ac) /agri, manu, serv/;
set f(ac) /lab, cap/;
alias(ac, acp), (a, ap), (f, fp);
table sam(ac, acp)
         agri   manu   serv    lab   cap    hh   total
agri    260     320    150                  635   1365
manu   45     390    390                  600    1725
serv    400     365    320                  385    1470
lab     200     250    400                          850
cap     160     400    210                          770
hh                             850   770            1620
total   1365    1725   1470    850    770    1620          ;
parameter rhoq(a)  /agri 0.2, manu 0.3,  serv 0.1/
          rhoVA(a)  /agri 0.25,  manu  0.3,  serv 0.1/;
parameter scaleAq(a)
          deltaq(a)
          scaleAVA(a)
          deltaVA(a)
          ia(a,ap)
          shrh(a)
          PA0(a)
QA0(a)
PVA0(a)
QVA0(a)
PINTA0(a)
QINTA0(a)
QINT0(a,ap)
QLD0(a)
QKD0(a)
QLS0
QKS0
WL0
WK0
YH0
QH0(a);

PA0(a)=1;
PVA0(a)=1;
PINTA0(a)=1;
WK0=1;
WL0=2;
QA0(a)=sam('total',a)/PA0(a);
QVA0(a)=sam('lab',a)/WL0+sam('cap',a)/WK0;
QINT0(a,ap)=sam(a,ap)/PA0(a);
QINTA0(a)=sum(ap,QINT0(ap,a));
ia(a,ap)=QINT0(a,ap)/QINTA0(ap);
QLS0=sum(a,sam('lab',a))/WL0;
QKS0=sum(a,sam('cap',a))/WK0;
QLD0(a)=sam('lab',a)/WL0;
QKD0(a)=sam('cap',a)/WK0;
deltaq(a)=PVA0(a)*QVA0(a)**(1-rhoq(a))/(PVA0(a)*QVA0(a)**(1-rhoq(a))+PINTA0(a)*QINTA0(a)**(1-rhoq(a)));
scaleAq(a)=QA0(a)/(deltaq(a)*QVA0(a)**rhoq(a)+(1-deltaq(a))*QINTA0(a)**rhoq(a))**(1/rhoq(a)) ;
deltaVA(a)=WL0*QLD0(a)**(1-rhoVA(a))/(WL0*QLD0(a)**(1-rhoVA(a))+WK0*QKD0(a)**(1-rhoVA(a)));
scaleAVA(a)=QVA0(a)/((deltaVA(a)*QLD0(a)**rhoVA(a))+((1-deltaVA(a))*QKD0(a)**rhoVA(a)))**(1/rhoVA(a));
YH0=WL0*QLS0+WK0*QKS0;
QH0(a)=sam(a,'hh')/PA0(a);
shrh(a)=(PA0(a)*QH0(a))/sum(ap,sam(ap,'hh'));
display PA0, PVA0, PINTA0, QA0, QVA0, QINTA0, QINT0, rhoq, rhoVA, scaleAq, deltaq,  deltaVA, scaleAVA, ia, shrh, QLD0, QKD0, QLS0, QKS0, WL0, WK0, YH0, QH0;
variable PA(a), PVA(a), PINTA(a), WL, WK, QA(a), QVA(a), QINTA(a), QINT(a,ap), QLD(a), QKD(a)
QLS, QKS, YH, QH(a);
equation EA(a), EB(a), EC(a), ED(a), EE(a), EF(a), EG(a,ap), EH(a), EI, EJ(a), EK(a),EL,EM;
EA(a)..  QA(a)=e= scaleAq(a)*(deltaq(a)*QVA(a)**rhoq(a)+(1-deltaq(a))*QINTA(a)**rhoq(a))**(1/rhoq(a));
EB(a)..  PVA(a)/PINTA(a) =e= (deltaq(a)/(1-deltaq(a)))*(QINTA(a)/QVA(a))**(1-rhoq(a));
EC(a)..   PA(a)*QA(a) =e= PVA(a)*QVA(a)+PINTA(a)*QINTA(a);
ED(a)..   QVA(a) =e= scaleAVA(a)*(deltaVA(a)*QLD(a)**rhoVA(a)+(1-deltaVA(a))*QKD(a)**rhoVA(a))** (1/rhoVA(a));
EE(a)..   WL/WK =e= (deltaVA(a)/(1-deltaVA(a)))*(QKD(a)/QLD(a))**(1-rhoVA(a));
EF(a)..   PVA(a)*QVA(a) =e= WL*QLD(a)+WK*QKD(a);
EG(a,ap)..   QINT(a,ap) =e= ia(a,ap)*QINTA(ap);
EH(ap)..  PINTA(ap) =e= sum(a, ia(a,ap)*PA(a));
EI..   YH =e= WL*QLS+WK*QKS;
EJ(a)..  PA(a)*QH(a)=e=shrh(a)*YH;
EK(a)..  QA(a)=e= sum(ap, QINT(a,ap)) +QH(a);
EL..  sum(a, QLD(a)) =e= QLS;
EM..  sum(a,QKD(a))=e= QKS;


PA.l(a)=PA0(a);
PVA.l(a)=PVA0(a);
PINTA.l(a)=PINTA0(a);
QA.l(a)=QA0(a);
QVA.l(a)=QVA0(a);
QINTA.l(a)=QINTA0(a);
QINT.l(a,ap)=QINT0(a,ap);
QLD.l(a)=QLD0(a);
QKD.l(a)=QKD0(a);
WK.l=1;
WL.l=1;
QLS.fx=QLS0;
QKS.fx=QKS0;
YH.l=YH0;
QH.l(a)=QH0(a);

model cge /all/;
solve cge using mcp;
QLS.fx=3*QLS0;
QKS.fx=9*QKS0;
model cge1 /all/;
solve cge1 using mcp;