How to minimise a value which is found by solving transcedental equations (using fsolve)?

I want to minimise the coefficient of variation of the y coordinate of the rocker chassis joint of the rocker bogie suspension system.

That is, the y coordinate of P5 in the figure.

The ground function is f(x,phi)=1+sin(x-phi)

P1 is at (0,f(0,phi))

phi is incremented to depict the motion of the rover as it moves along the ground. (initially it is 0)

Each of the points P2,P3,P4 and P5 are found sequentially using fsolve. (I am not able to express P5 as a function of P1 analytically. This is why I don't know how to minimise it.)

The equations are obtained from coordinate geometry as the triangles made by P1,P3,P2 and P3,P5,P4 are fully determined. It is important to note that the angle between P2,P3,P5 is not constant.

Once the y coordinate of P5 is found, I will store it in an array and loop it over phi. Then the coefficient of variation of all the y values that have been collected is computed.

After this comes the minimisation part that I'm stuck at.

The code for phi=0 is attached below. The steps after this (looping phi and storing values) are elementary until the minimisation.

Note that the code uses the name X instead of P.

%% Function to find points of a rocker bogie traversing a function given the link lengths and angles

function y=x5solver(L,ground)

%% Defining initial stuff
%ground=@(x,phi)1+sin(x-phi);
X(1,1)=0;
d12=sqrt(L(1)^2+L(2)^2-2*L(1)*L(2)*cos(L(5)));
d34=sqrt(L(3)^2+L(4)^2-2*L(3)*L(4)*cos(L(6)));
d13=L(1);
d23=L(2);
d35=L(3);
d45=L(4);


options=optimset('Display','off');

%% Phase = 0, Solving for X(2,:)
phase=0;
f=@(x)ground(x,phase);
X(1,2)=f(X(1,1));

fun=@(x)ground_distance_solver(x,X(1,:),d12,f);
x0=[d12 0];
X(2,:)=fsolve(fun,x0,options);

%% Solving for X(3,:)
fun=@(x)two_distance_solver(x,X(1,:),X(2,:),d13,d23);
x0=[(X(1,1)+X(2,1))/2 (X(1,2)+X(2,2))/2+3];
X(3,:)=fsolve(fun,x0,options);

%% Solving for X(4,:)
fun=@(x)ground_distance_solver(x,X(3,:),d34,f);
x0=[d34 0];
X(4,:)=fsolve(fun,x0,options);

%% Solving for X(5,:)
fun=@(x)two_distance_solver(x,X(3,:),X(4,:),d35,d45);
x0=[(X(3,1)+X(4,1))/2 (X(3,2)+X(4,2))/2+3];
X(5,:)=fsolve(fun,x0,options);

%%
y=X;
%% Display for testing
% disp(X);
% xplot=linspace(-1,100,1000);
% yplot=f(xplot);
% plot(xplot,yplot);
% axis([-1 100 -1 100]);
% hold on 
% scatter(X(:,1),X(:,2));
% hold off

Thanks for your time and effort.