stimulate small strain

i am trying to stimulate the small strain level, i set the loadrate, timestep,and max_vel very small. but i cannot get a constant small strain stiffness, the curve of stress-strain always has fluctuation. i cannot figure it out . here is my code

FrictionAngle = 35

Damp = 0.25

#Confining variables
ConfPress = -2.0e5

#Loading control
LoadRate = -0.001
O.load('calm.yade.bz2')

########################################
#import necessary packages
from yade import pack,plot,os,timing

triax=O.engines[4]
#################################
#####Defining triaxil engines####
#################################

triax=TriaxialStressController(
    thickness = 0.001,
    internalCompaction = False, # If true the confining pressure is generated by growing particles
    max_vel = 1e-20,
    stressMask = 5,
    goal1 = ConfPress,
    goal2 = LoadRate,
    goal3 = ConfPress,

)

ini_e22a = -triax.strain[1]
ini_e22b = -triax.strain[1]

#O.usesTimeStepper=True

O.dt=1e-6

newton=NewtonIntegrator(damping=Damp)

setContactFriction(radians(FrictionAngle))

O.trackEnergy=True

O.timingEnabled=True

###engine
O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    InteractionLoop(
     [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
     [Ip2_FrictMat_FrictMat_MindlinPhys()],
     [Law2_ScGeom_MindlinPhys_Mindlin()]
    ),
    #GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    triax,
    newton,
    PyRunner(iterPeriod=1,command='endCheck()',label='check'),
    PyRunner(command='addPlotData()',iterPeriod=10,label='record'),
    #VTKRecorder(iterPeriod=100,recorders=['all'],fileName='./post/p1-'),
    TriaxialStateRecorder(iterPeriod=100,file='WallStresses3')
]
# Simulation stop conditions defination
def endCheck():
    unb=unbalancedForce()
    e22=-triax.strain[1]
    if abs(e22)>1.5e-5:
       O.pause()
       O.save('e0.05.yade.bz2')
# collect history of data
def addPlotData():

    global ini_e22a
    global ini_e22b

    e22=-triax.strain[1]
    unb = unbalancedForce()
    mStress = -(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3.0
    volume = (O.bodies[1].state.pos[0]-O.bodies[0].state.pos[0])*(O.bodies[3].state.pos[1]-O.bodies[2].state.pos[1])*(O.bodies[5].state.pos[2]-O.bodies[4].state.pos[2])
    Porosity=1-sum(4*pi*b.shape.radius**3/3 for b in O.bodies if isinstance(b.shape,Sphere))/volume

    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],
                 ub=unbalancedForce(),
                 dstress=(-triax.stress(triax.wall_top_id)[1])-(-triax.stress(triax.wall_right_id)[0]),
                 disx=triax.width,
                 disy=triax.height,
                 disz=triax.depth,
   i=O.iter,
                 porosity=Porosity,
                 cnumber=avgNumInteractions(),
                 ke=utils.kineticEnergy(),
                 totale=O.energy.total()
                 )

    print ('unbalanced force: %f, mean stress: %f, s11: %f, s22: %f, s33: %f, coordination number: %f, porosity: %f' %(unb,mStress,-triax.stress(triax.wall_right_id)[0],-triax.stress(triax.wall_top_id)[1],-triax.stress(triax.wall_front_id)[2],avgNumInteractions(),Porosity))
    plot.saveDataTxt('loadinglog3.txt.bz2')
##################################################
######Defining output for postprocessing##########
##################################################

plot.plots={'e22':('s22',None,'ev')}

plot.plot()