function Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt)

Mat_G=[
    sparse(PD(Loadi))-mPD(Loadi)-bigM*PDbi-uPD(Loadi);
    sparse(PD(Loadi))-mPD(Loadi)+bigM*PDbi+lPD(Loadi);
    sparse(QD(Loadi))-mQD(Loadi)-bigM*QDbi-uQD(Loadi);
    sparse(QD(Loadi))-mQD(Loadi)+bigM*QDbi+lQD(Loadi);
    Volt'-mVolt'-bigM*Vbi-uVolt;
    Volt'-mVolt'+bigM*Vbi+lVolt;
    Vbi;
    PDbi;
    QDbi;
    (0.1-Vbi).*Vbi;
    (0.1-PDbi).*PDbi;
    (0.1-QDbi).*QDbi;
    ];
end