Thank you for your comment.
I agree that return np.array([g0]) is suspicious, but to my surprise, it does not.
For instance if I have:
g0 = 1 / con - beta * alpha

instead of
g0 = 1 / con - beta * alpha * torch.sum(omega5 * (
            1 / con_plus * aplus * x[0]**(alpha-1) * ls_plus**(1-alpha)))

In this case, PyTorch correctly computes Jacobian. I meant the Jacobian is numerically equivalent to the one from the IPOPT’s derivative checker.
This is why I am confused…
How do you see the issue?