def objFunc(p):
    x0Val, z0Val, x1Val = p
    C = (1,1,1)
    A = (x0Val, 0, z0Val)
    B = (x1Val, 0, 0)
    #result already divided by 2
    return 4*((1-x0Val)**2+1+(1-z0Val)**2)**0.5 + 2*((x0Val-x1Val)**2+z0Val**2)**0.5 + x1Val


sol = minimize_constrained(objFunc, [(0, 1)]*3, (0.7, 0.5, 0.6))
print ("x0, z0, x1 = ", sol)
print ("total length = ", objFunc(sol))


def getGraphic(x0Val, z0Val, x1Val):
    #[(1, y, z) for y,z in itertools.product(*[[-1, 1]]*2)]
    C1 = (1,1,1)
    C2 = (1,-1,1)
    C3 = (1,1,-1)
    C4 = (1,-1,-1)
    C5 = (-1,1,1)
    C6 = (-1,-1,1)
    C7 = (-1,1,-1)
    C8 = (-1,-1,-1)
    
    A1 = (x0Val, 0, z0Val)
    A2 = (x0Val, 0, -z0Val)
    B = (x1Val, 0, 0)
    B2 = (-x1Val, 0, 0)
    A3 = (-x0Val, 0, z0Val)
    A4 = (-x0Val, 0, -z0Val)
    
    allLines = (
        [(C1, A1), (C2, A1), (C3, A2), (C4, A2), (A1, B), (A2, B),
         (B, B2),
        (C5, A3), (C6, A3), (C7, A4), (C8, A4), (A3, B2), (A4, B2)]
    )
    
    g = Graphics()
    g += points([C1,C2,C3,C4,C5,C6,C7,C8], size=20, color='blue')
    for t in allLines: g += line(t, color='green', radius=0.015)
    
    return g

getGraphic(*sol).show()

#Output:
#x0, z0, x1 =  (0.7114348459256075, 0.5006602988247546, 0.4225927055984936)
#total length =  6.196153088395327