The answer to the 2 questions:
1.
What would be the mathematical way to obtain the transfer function for this system and find gains from it?
let reference signal = r;
position(rad) = y;
then the transfer function from r to y (I calculated it algebraically):
Let parallel PI form: PIp = Kpp + Kip/s, series PI form PIs = Kps(1+Kis/s)
in discrete: ideal(series) PI = P*(1+Iα(z))
parallel PI:
(y/r)(s) = ((Kpp*s^2) + Kip )/ (Kip + (1-Kpp)*s^2)
series (aka ideal) PI:
(y/r)(s) = ((Kps*Kis*s^2) + Kis)/ (Kis + (1-Kps*Kis)*s^2)
I can do it in discrete time / frequency case ( in z-domain) as well if you accept the my answer and hire contact me for.
by Tustin/Bilinear continuous to discrete: s -> (2/Ts)*(z-1)/(z+1)
Then the gains are Kpp, Kip, Kps, Kis you can find in most simpler and clear way from the transfer function desired behaviour by zeroes poles placement or via settling time, rise time, margins, overshoot, undershoot formulas for the resultant second order closed loop system. (let me know if you need to construct the exact expressions/formulas further).
2.
How would one conduct a stability analysis for this system using matlab transfer function and bode plots when the gains to the PI controller are unknown?
figure
w = logspace(2,6,1000);
bode(closed_loop_tr_f_in_2_vectors_form,w),grid
In TCE NCE MPPL Matlab/Octave you can check whether the resultant closed loop system is stable (and also is stabilizable in Octave already, not yet in Matlab) by just invoking built-in function:
In TCE NCE MPPL transfer function is presented as 2 vectors b for numerator and a for denumerator coefficients.
%for instance for parallel PI form: (for seria(ideal) the way is the same, just the transfer function denominator and numerator coefficient to substitute accordingly):
isstable(tf2ss([Kpp 0 Kip] [1-Kpp 0 Kip]))
With Bode:
closed_loop_tr_f_in_2_vectors_form = tf([Kpp 0 Kip] [1-Kpp 0 Kip])
figure
w = logspace(2,6,1000);
bode(closed_loop_tr_f_in_2_vectors_form,w),grid
%In Octave also isstabilizable(y/r(s))
isstabilizable(tf2ss([Kpp 0 Kip] [1-Kpp 0 Kip]))
Another way is to find and check the zeros and poles locations and so stability from the transfer function mathematically or via root locus analysis (can do it if you hire me).