How to reach target confining stress in triaxial test

Hi,
I am using [1] for simulating a triaxial test. The following MWE is from [1] and only some parameters are modified (damp, targetPorosity, stabilityThreshold).
My question is, after the confining process is completed, the confining pressure doesn't reach our desired pressure which is 10kPa, the mean stress and other 3 stresses are:
meanStress: 9997.34994518; top: 9649.20616727 ;right: 9842.19217112 ;front :10510.4130631.
It is very nice for meanStress, but the value of other 3 directions are far from the target goal, and this results in the stress-strain curve doesn't start from origin point.
Here is the MWE, it takes about 30 seconds to run it.
#####
from yade import pack, plot

nRead=readParamsFromTable(
 num_spheres=2000,
 compFricDegree = 30,
 key='_triax_base_',
 unknownOk=True
)
from yade.params import table

num_spheres=table.num_spheres
key=table.key
targetPorosity = 0.7
compFricDegree = table.compFricDegree
finalFricDegree = 30
rate=-0.01
damp=0.4
stabilityThreshold=0.001
young=5e7

confinement=10e3

mn,mx=Vector3(0,0,0),Vector3(0.7,1.4,0.7)

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))

## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

Gl1_Sphere.quality=3

triax=TriaxialStressController(

 maxMultiplier=1.+4e-3,
 finalMaxMultiplier=1.+4e-4,
 thickness = 0,

 stressMask = 7,
 internalCompaction=True,
)
newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
 newton
]

Gl1_Sphere.stripes=0
if nRead==0: yade.qt.Controller(), yade.qt.View()

triax.goal1=triax.goal2=triax.goal3=-confinement

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print 'unbF:',unb,' meanStress: ',-triax.meanStress,'top:',-triax.stress(triax.wall_top_id)[1]
  if unb<stabilityThreshold and abs(-confinement-triax.meanStress)/confinement<0.001:
 break

print "### Isotropic state saved ###"
##### below is what I tried, but it seems doesn't work
# while 1:
# O.run(100, True)
# print ' meanStress: ',-triax.meanStress,'top:',-triax.stress(triax.wall_top_id)[1],'front:',-triax.stress(triax.wall_front_id)[2],\
# 'right:',-triax.stress(triax.wall_right_id)[0]
# if abs(-confinement-triax.meanStress)/confinement<0.001 and \
# abs(-confinement-triax.stress(triax.wall_right_id)[0])/confinement<0.01 and \
# abs(-confinement-triax.stress(triax.wall_top_id)[1])/confinement<0.01 and \
# abs(-confinement-triax.stress(triax.wall_front_id)[2])/confinement<0.01:
# break

import sys
while triax.porosity>targetPorosity:
 # we decrease friction value and apply it to all the bodies and contacts
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 # while we run steps, triax will tend to grow particles as the packing
 # keeps shrinking as a consequence of decreasing friction. Consequently
 # porosity will decrease
 O.run(500,1)

#O.save('compactedState'+key+'.yade.gz')
print "### Compacted state saved ###"

print 'meanStress: ',-triax.meanStress, 'top', -triax.stress(triax.wall_top_id)[1], 'right',-triax.stress(triax.wall_right_id)[0],'front',-triax.stress(triax.wall_front_id)[2]

triax.internalCompaction=False
setContactFriction(radians(finalFricDegree))
triax.stressMask = 5
triax.goal2=rate
triax.goal1=-confinement
triax.goal3=-confinement

newton.damping=0.4
O.saveTmp()

def history():
 plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
   ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
   s11=-triax.stress(triax.wall_right_id)[0],
   s22=-triax.stress(triax.wall_top_id)[1],
   s33=-triax.stress(triax.wall_front_id)[2],
   dev=-triax.stress(triax.wall_top_id)[1]-10000,
   i=O.iter)

def stopIfDamaged():
 if -triax.strain[1]>0.1:
  O.pause()
  plot.saveDataTxt('Frict_10kpa_num{}_CaAtGap'.format(num_spheres))

if 1:
 O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+[PyRunner(iterPeriod=20,command='stopIfDamaged()',label='checkDamage')]+O.engines[5:7]
plot.plots={'e22':('s11','s22','s33'),'e22 ':('dev')}
plot.labels={'e22':'$\epsilon_{a}$','s11':'$\sigma_{11}$','s22':'$\sigma_{22}$','s33':'$\sigma_{33}$'}
plot.plot()
###########
The commented part of code is what I tried to let the stress reach our target value. but it seems doesn't work since the top, front and right stress doesn't change even after run it for several hours.

Could you please give some instruction about how to reach target confining stress?
Thanks
Leonard
[1]https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py