• and are mass and inertia moment.
• Geometry. : hub to saddle; : wheel radius; : crank length.
• Normal forces: is the total normal force. is the force in saddle (and hub). is in the front pedal, in rear pedal. σ is the proportion of that is loaded on the saddle; β is the proportion of the remaining that is loaded in the front pedal.
• Tangential forces. is in front pedal, in rear pedal. is the tangential force on the wheel.
• are the horizontal and vertical displacements of the center of gravity.
• is the rotation of the unicycle with respect to the hub.
• θ is the rotation of the cranks.

In the horizontal axis thus

In the vertical axis so The normal forces are: and

I've used sage to manipulate these equations. This is the code, in case anybody wants to help:

# Mass, g and inertia moment
m, g, I = var('m, g, I')

# Geometry. L: hub to saddle; R: wheel radius; r: crank length
L, R, r = var('L, R, r')

# Normal forces.  S is in saddle (and hub)). N1 is in front pedal, N2
# in rear pedal.  Sigma is the proportion of N that is loaded on the
# saddle; beta is the proportion of the remaining N that is loaded in
# the front pedal.
N, S, N1, N2, sigma, beta = var('N, S, N1, N2, sigma, beta')

# Tangential forces.  T1 is in front pedal, T2 in rear pedal. 
T, T1, T2 = var('T, T1, T2')

t = var('t')
phi = function('phi', t)
theta = function('theta', t)

x(t) = theta*R + L*sin(phi)
ddx(t) = diff(x, t, 2)

y(t) = L - L*cos(phi)
ddy(t) = diff(y, t, 2)

S = sigma*N
N1 = beta*(1-sigma)*N
N2 = (1-beta)*(1-sigma)*N

eqT1 = solve((T1*(R + r*cos(theta)) + T2*(R + r*cos(theta+pi)) + T*R
              -N1*r*sin(theta) - N2*r*sin(theta + pi) == 0), T)[0]
eqN = solve((m*(g + ddy(t)) == S + N1 + N2), N)[0]

eqT = (m*ddx(t) == T1 + T2 + T)(T = eqT1.right())
eqT = eqT(N = eqN.right())

eqI = (-I*diff(phi, 2) ==
        -S*L*sin(phi) + T*L*cos(phi)
        -N1*(L*sin(phi) - r*sin(theta))
        -N2*(L*sin(phi) - r*sin(theta + pi))
        +T1*(L*cos(phi) - r*cos(theta))
        +T2*(L*cos(phi) - r*cos(theta + pi)))(T = eqT1.right())
eqI = eqI(N = eqN.right())

# Substitute T1=T2=0 in eqT and eqI, and then edit to replace cos(phi)
# with 1 and sin(phi) with phi. 

eqISimp = diff(phi(t), t, 2)*I == ((1 - 2*beta)*r - (1 - 2*beta)*r*sigma)*sin(theta(t))*L*((m*phi(t)*diff(phi(t), t, 2) + m*diff(phi(t), t, 1)^2)*L + g*m)/R - (beta - 1)*(1 - sigma)*(phi(t)*L + r*sin(theta(t)))*((m*phi(t)*diff(phi(t), t, 2) + m*diff(phi(t), t, 1)^2)*L + g*m) + beta*(1 - sigma)*(phi(t)*L - r*sin(theta(t)))*((m*phi(t)*diff(phi(t), t, 2) + m*diff(phi(t), t, 1)^2)*L + g*m) + sigma*phi(t)*L*((m*phi(t)*diff(phi(t), t, 2) + m*diff(phi(t), t, 1)^2)*L + g*m)

eqTSimp = m*(diff(theta(t), t, 2)*R + diff(phi(t), t, 2)*L - phi(t)*diff(phi(t), t, 1)^2*L) == ((1 - 2*beta)*r*sigma - (1 - 2*beta)*r)*sin(theta(t))*((m*phi(t)*diff(phi(t), t, 2) + m*diff(phi(t), t, 1)^2)*L + g*m)/R

#eqISimp = (-diff(phi(t), t, 2)*I == ((beta*r - beta*r*sigma)*sin(theta(t)) - ((1 - beta)*r*sigma + (beta - 1)*r)*cos(theta(t)))*L*((m*phi(t)*diff(phi(t), t, 2) + m*diff(phi(t), t, 1)^2)*L + g*m)/R - beta*(1 - sigma)*(phi(t)*L - r*sin(theta(t)))*((m*phi(t)*diff(phi(t), t, 2) + m*diff(phi(t), t, 1)^2)*L + g*m) + (beta - 1)*(1 - sigma)*(phi(t)*L - r*cos(theta(t)))*((m*phi(t)*diff(phi(t), t, 2) + m*diff(phi(t), t, 1)^2)*L + g*m) - sigma*phi(t)*L*((m*phi(t)*diff(phi(t), t, 2) + m*diff(phi(t), t, 1)^2)*L + g*m))
#eqTSimp = m*(diff(theta(t), t, 2)*R + diff(phi(t), t, 2)*L - phi(t)*diff(phi(t), t, 1)^2*L) == ((beta*r - beta*r*sigma)*sin(theta(t)) - ((1 - beta)*r*sigma + (beta - 1)*r)*cos(theta(t)))*((m*phi(t)*diff(phi(t), t, 2) + m*diff(phi(t), t, 1)^2)*L + g*m)/R