Minimize cost function (4-d integral)

I want to integrate this cost function symbolically in order to later find its minimum using fmincon():

where K:

and , depend on , , , , , , θ. However, the int() function fails to integrate. The example 4-D Integral of Sphere is not suitable, since there the integral is calculated numerically.

[Dyuman Joshi]

@Arina Pirogova, The images are not of clear resolution/quality than the ones I edited to add in the question description.

Add them as images rather than files.

[Arina Pirogova]

I did, but they are not shown in preview.

[Sam Chak]

@Arina Pirogova, To attach an image, you have to click this icon from the INSERT submenu.

[Sam Chak]

@Arina Pirogova, Usually, we use fmincon numerically.

fun = @(x) (x - 1).^2;
x0  = 0;
x   = fmincon(fun, x0, [], [])

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

x = 1.0000

But if your cost function is described symbolically, this happens:

syms x
f(x) = (x - 1)^2

f(x) = 

x0  = 0;
x   = fmincon(f, x0, [], [])

Error using fmincon
FMINCON requires all values returned by functions to be of data type double.

[Arina Pirogova]

@Sam Chak of course, I have for x_0 some values of psi_i, phi_i i=1,2. But I don’t understand how to represent a cost function in order to minimize it.

[Torsten]

But I don’t understand how to represent a cost function in order to minimize it.

Either by defining a function handle or a function with p = [psi1,psi2,phi1,phi2] as input vector from which you deduce the value of the cost function in this specific point.

Look at the examples provided for "fmincon", e.g.

[Arina Pirogova]

Do you mean like that? Error: Not enough input arguments.

    syms t theta ddtheta ddx ddz psi1 psi2 phi1 phi2
 A=[1 cos(phi1) cos(phi2);0 sin(phi1) sin(phi2);...
    0 sin(psi1)/abs(sin(psi1+phi1)) sin(psi2)/abs(sin(psi2+phi2))];
T_theta=[cos(theta) sin(theta) 0;-sin(theta) cos(theta) 0;0 0 1];
B=[ddx; ddz+k*g; mu*ddtheta];
U=inv(A)*T_theta*B;
u0=abs(U(1));
u1=abs(U(2));
u2=abs(U(3));
UU=simplify(u0+u1+u2);
func=matlabFunction(UU,'Vars',[theta,ddtheta,ddx,ddz,psi1,psi2,phi1,phi2]);
I3=@(theta) integral3(@(ddtheta,ddx,ddz) func(theta,ddtheta,ddx,ddz,psi1,psi2,phi1,phi2),-delta/mu,delta/mu,-delta,delta,-delta,delta);
f = @(psi1,psi2,phi1,phi2) integral(I3,-pi,pi);
psi1_0 = 3*pi/4;
psi2_0 = 0;
phi1_0 = 5*pi/12;
phi2_0 = 7*pi/4;
x_0 = [psi1_0,psi2_0,phi1_0,phi2_0]; 
lb = [0,pi,0,pi];
ub =[pi,2*pi,pi,2*pi];
fmincon(f,x_0,[],[],[],[],lb,ub,@non_linear);
%Nonlinear constraints:
function [c,ceq]=non_linear(x)
  c = [-abs(tan(x(1)+x(3))+0.2 
-abs(tan(x(2)+x(4))+0.2];
  ceq =[];
end

[Sam Chak]

@Arina Pirogova, I tried to fix what I can. If possible, put some annotations or comments.

%% Symbolic work section
syms t theta ddtheta ddx ddz psi1 psi2 phi1 phi2 g k mu
%                                                ^ ^ ^  missing variables

% Some matrix A
A       = [1 cos(phi1) cos(phi2);
           0 sin(phi1) sin(phi2); 
           0 sin(psi1)/abs(sin(psi1+phi1)) sin(psi2)/abs(sin(psi2+phi2))];

% Some matrix theta (looks like a Rotation Matrix)
T_theta = [cos(theta) sin(theta) 0;
          -sin(theta) cos(theta) 0;
                    0          0 1];

% Some matrix B
B       = [ddx; 
           ddz+k*g; 
           mu*ddtheta];

% Maybe some state equations?
U       = inv(A)*T_theta*B;
u0      = abs(U(1));
u1      = abs(U(2));
u2      = abs(U(3));
UU      = simplify(u0 + u1 + u2);

%% Numerical work section
func    = matlabFunction(UU, 'Vars', [theta, ddtheta, ddx, ddz, psi1, psi2, phi1, phi2, g, k, mu]);

% don't know what 'delta' is, but the anon function takes it.
I3      = @(theta) integral3(@(ddtheta, ddx, ddz) func(theta, ddtheta, ddx, ddz, psi1, psi2, phi1, phi2), -delta/mu, delta/mu, -delta, delta, -delta, delta)

I3 = function_handle with value:
    @(theta)integral3(@(ddtheta,ddx,ddz)func(theta,ddtheta,ddx,ddz,psi1,psi2,phi1,phi2),-delta/mu,delta/mu,-delta,delta,-delta,delta)

% The cost function
f       = @(psi1, psi2, phi1, phi2) integral(I3, -pi, pi);

psi1_0  = 3*pi/4;                           % for x_0
psi2_0  = 0;                                % for x_0
phi1_0  = 5*pi/12;                          % for x_0
phi2_0  = 7*pi/4;                           % for x_0
x_0     = [psi1_0,psi2_0,phi1_0,phi2_0];    % initial guess

lb      = [0,pi,0,pi];                      % lower bound
ub      = [pi,2*pi,pi,2*pi];                % upper bound

%% call fmincon solver
[x, fval] = fmincon(f, x_0, [], [], [], [], lb, ub, @non_linear)

delta requires RF PCB Toolbox.

Error in solution>@(theta)integral3(@(ddtheta,ddx,ddz)func(theta,ddtheta,ddx,ddz,psi1,psi2,phi1,phi2),-delta/mu,delta/mu,-delta,delta,-delta,delta) (line 28)
I3 = @(theta) integral3(@(ddtheta, ddx, ddz) func(theta, ddtheta, ddx, ddz, psi1, psi2, phi1, phi2), -delta/mu, delta/mu, -delta, delta, -delta, delta)

Error in integralCalc/iterateScalarValued (line 314)
fx = FUN(t);

Error in integralCalc/vadapt (line 132)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 87)
Q = integralCalc(fun,a,b,opstruct);

Error in solution>@(psi1,psi2,phi1,phi2)integral(I3,-pi,pi) (line 30)
f = @(psi1, psi2, phi1, phi2) integral(I3, -pi, pi); % <-- cost function

Error in fmincon (line 563)
initVals.f = feval(funfcn{3},X,varargin{:});

%% Nonlinear constraints:
function [c, ceq] = non_linear(x)
    c   = [-abs(tan(x(1) + x(3)) + 0.2);    % <-- unsure if the bracket is placed correctly
           -abs(tan(x(2) + x(4)) + 0.2)];   % <-- unsure if the bracket is placed correctly
    ceq = [];
end

[Arina Pirogova]

%Constants:

 mu = 55;
    delta = 1;
    k = 0.142857;
    g = 9.8;

if you assign these values to variables, then fmincon throws an error: Input function must return ‘double’ or ‘single’ values. Found 'sym'.

I tried to substitute x_0 into the function f to check its functionality. It didn't work for the same reason. So, please explain how to pass the value psi_i, phi_i to the function f, since there is only I3, and their dependence is not specified.

[Torsten]

Check whether this is really what you want. The construction is quite complicated, and the evaluation of your function F takes a while.

%% Symbolic work section
syms theta ddtheta ddx ddz psi1 psi2 phi1 phi2 
%Constants:
mu = 55;
delta = 1;
k = 0.142857;
g = 9.8;

% Some matrix A
A       = [1 cos(phi1) cos(phi2);
           0 sin(phi1) sin(phi2); 
           0 sin(psi1)/abs(sin(psi1+phi1)) sin(psi2)/abs(sin(psi2+phi2))];

% Some matrix theta (looks like a Rotation Matrix)
T_theta = [cos(theta) sin(theta) 0;
          -sin(theta) cos(theta) 0;
                    0          0 1];

% Some matrix B
B       = [ddx; 
           ddz+k*g; 
           mu*ddtheta];

% Maybe some state equations?
U       = inv(A)*T_theta*B;
u0      = abs(U(1));
u1      = abs(U(2));
u2      = abs(U(3));
UU      = simplify(u0 + u1 + u2);

%% Numerical work section
func    = matlabFunction(UU, 'Vars', [theta, ddtheta, ddx, ddz, psi1, psi2, phi1, phi2]);

I3      = @(theta,psi1,psi2,phi1,phi2) integral3(@(ddtheta, ddx, ddz) func(theta, ddtheta, ddx, ddz, psi1, psi2, phi1, phi2), -delta/mu, delta/mu, -delta, delta, -delta, delta);
% The cost function
f       = @(psi1, psi2, phi1, phi2) integral(@(theta)I3(theta,psi1,psi2,phi1,phi2), -pi, pi, 'ArrayValued',1);
F       = @(x)f(x(1),x(2),x(3),x(4));
psi1_0  = 3*pi/4;                           % for x_0
psi2_0  = 0;                                % for x_0
phi1_0  = 5*pi/12;                          % for x_0
phi2_0  = 7*pi/4;                           % for x_0
x_0     = [psi1_0,psi2_0,phi1_0,phi2_0];    % initial guess
lb      = [0,pi,0,pi];                      % lower bound
ub      = [pi,2*pi,pi,2*pi];                % upper bound

%% call fmincon solver
[x, fval] = fmincon(F, x_0, [], [], [], [], lb, ub, @non_linear)
%% Nonlinear constraints:
function [c, ceq] = non_linear(x)
    c   = [-abs(tan(x(1) + x(3))) + 0.2;    % <-- unsure if the bracket is placed correctly
           -abs(tan(x(2) + x(4))) + 0.2];   % <-- unsure if the bracket is placed correctly
    ceq = [];
end