How to reduce unb while keeping state constant

Hi,
I read from the book <Particulate Discrete Element Modelling> in which the author suggests that in servo-controlled simulations once the desired specific stress state is attained, it is good practice to halt all boundary motion and run through a number of DEM calculation cycles to ensure that the system is in equilibrium at that required stress level.

So I want to halt the 6 walls in the triaxial compression test and then run some steps. Here I tied in 2 ways:

The first approach I came up with was fixing all the walls (please see script2-approach 1)

After I O.run(10000), I found that: 1)the unbalancedForce decreased which is expected, 2) triax.porosity is the same which is good. 3) However, the stress on each wall seems to become Vector3(0,0,0), which means the stress state of the sampel is not remained.

The second approach I tried was to keep the stress of all the walls constant as confining pressure (please see script2-approach 2).
After I O.run(10000), I found that: 1) unb decreased which is expected, 2) triax.porosity changed which is not desired because it means the dense state of the sample changed. 3) the stress on each wall get more precisely close to the confining stress which is good.

Both the two ways couldn't keep the current state , do you have any ideas to let simulations run but keep current state?
---
Here is the MWE (modified from Bruno's script) if you want to reproduce this. The MWE is in two-steps: 1st for making sample and 2nd for reloading sample and checking.
######script1#######
from __future__ import division
from yade import pack, plot
num_spheres=7000#
targetPorosity = 0.44
compFricDegree = 30
finalFricDegree = 18
rate=-0.01
damp=0.6
stabilityThreshold=0.001
young=2e8
confinement=100e3

mn,mx=Vector3(0,0,0),Vector3(0.07,0.14,0.07)

MatWall=O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
MatSand = O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.3,num_spheres,False, 0.95,seed=1)

O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
triax=TriaxialStressController(
    maxMultiplier=1.+3e7/young,
    finalMaxMultiplier=1.+16e4/young,
    thickness = 0,
    stressMask = 7,
    internalCompaction=True,
)
newton=NewtonIntegrator(damping=damp)

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    triax,
    TriaxialStateRecorder(iterPeriod=100,file='WallStresses'),
    newton,
]

Gl1_Sphere.stripes=0
triax.goal1=triax.goal2=triax.goal3=-confinement

while 1:
  O.run(1000, True)
  #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  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.0001:
      break

print "### Isotropic state saved ###"

import sys
while triax.porosity>targetPorosity:
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 O.run(500,1)

O.save('sample.yade.gz')
print "### Compacted state saved ###"

######script2#######
from yade import pack, plot

O.load('sample.yade.gz')
confinement=100e3
triax=TriaxialStressController(
    thickness = 0,
    stressMask = 5,
    internalCompaction=False,
)
newton=NewtonIntegrator(damping=0.4)

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True),Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    triax,
    TriaxialStateRecorder(iterPeriod=1000,file='WallStresses'),
    newton,
]
triax.stressMask=7
O.step()

## Try approach 1
triax.wall_top_activated=False
triax.wall_bottom_activated=False
triax.wall_front_activated=False
triax.wall_back_activated=False
triax.wall_left_activated=False
triax.wall_right_activated=False

## Try approach 2
# triax.goal1=triax.goal2=triax.goal3=-confinement

print "before run:",triax.stress(triax.wall_front_id),triax.stress(triax.wall_top_id),triax.stress(triax.wall_left_id)
print "porosity:",triax.porosity,"unb:",unbalancedForce()
O.run(10000,True)
print "after run:",triax.stress(triax.wall_front_id),triax.stress(triax.wall_top_id),triax.stress(triax.wall_left_id)
print "porosity:",triax.porosity,"unb:",unbalancedForce()

Thanks
Leonard