Thank you for your guidance. To avoid confusion with the variables in the original question, I replaced z with δ and set distinct upper bounds for x and y (u_x and u_y). Here’s my adjusted code:
model = Model()
@variable(model, 0 <= x <= u_x)
@variable(model, 0 <= y <= u_y)
@variable(model, δ, Bin) # Using delta as the binary variable
# δ == 0 --> 0 <= x <= 0
# δ == 1 --> 0 <= x <= u_x
@constraint(model, x <= u_x * δ)
# δ == 0 --> 0 <= y <= u_y
# δ == 1 --> u_y <= y <= u_y
@constraint(model, u_y - y <= u_y * (1 - δ))
If there’s anything inconsistent with your approach, please let me know.