optimization expression includes an integration

I am trying an optimization problem in which the expression of the objective function includes an integral.

It is obvious that the \sigma equal to one results in the optimium solution. I want to use the optimization toolbox to get this result with an initial \sigma equal to, say, 10.

I wrote the following code.

g1 = @(x,c) (exp(-(0.5*(x./c).^2))./sqrt(2*pi*c^2));
c = optimvar("c",1,1,'Type','continuous','LowerBound',0.1,'UpperBound',10);
prob = optimproblem('Objective', (0.5 - integral(@(x)g1(x,c),0, 10)).^2);
[solf,fvalf,eflagf,outputf] = solve(prob)

The following error is generated.

Error using integralCalc>finalInputChecks (line 544)

Input function must return 'double' or 'single' values. Found

'optim.problemdef.OptimizationExpression'.

I have two questions:

1, Am I coding the problem properly/correctly?

2, If the code is basically correct, how can I solve the error?

Thank you.

Answer 1

[Torsten]

c = optimvar("c",1,1,'Type','continuous','LowerBound',0.1,'UpperBound',10);
g1 = @(x,c) exp(-0.5*(x./c).^2)./sqrt(2*pi*c^2);
obj = fcn2optimexpr(@(c)abs(0.5-integral(@(x)g1(x,c),0,10)).^0.5,c);
prob = optimproblem('Objective', obj);
show(prob)

  OptimizationProblem : 

	Solve for:
       c

	minimize :
       arg1

       where:

           anonymousFunction1 = @(c)abs(0.5-integral(@(x)g1(x,c),0,10)).^0.5;
           arg1 = anonymousFunction1(c);


	variable bounds:
       0.1 <= c <= 10

x0.c = 10;
options = optimoptions("fmincon",OptimalityTolerance=1e-18);
[solf,fvalf,eflagf,outputf] = solve(prob,x0,Options=options)

Solving problem using fmincon.

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.

solf = struct with fields:
    c: 1.0115

fvalf = 0

eflagf = 
    OptimalSolution

outputf = struct with fields:
              iterations: 6
               funcCount: 35
         constrviolation: 0
                stepsize: 6.9192e-06
               algorithm: 'interior-point'
           firstorderopt: 0
            cgiterations: 18
                 message: 'Local minimum found that satisfies the constraints....'
            bestfeasible: [1x1 struct]
     objectivederivative: "finite-differences"
    constraintderivative: "closed-form"
                  solver: 'fmincon'

%The function showing the error is very flat - thus c = 1 is unlikely as
%result.
y = 0.01:0.01:3;
G1 = arrayfun(@(y)abs(0.5-integral(@(x)g1(x,y),0,10)).^0.5,y);
plot(y,G1)
grid on

[Sukvasant Tantikovt]

Hi Torsten,

Thank you for your help.

The solution seems to be very sensitivty to the upper bound value. If this value is changed to 30, the code produces c equal to 1.5 747. If the upper bound is 20, c = 0.4068 is obtained. Do you have any idea?

[Torsten]

I plotted the error function for your problem. It's almost 0 for c in the interval [0 3] and too flat to get good convergence. Most probably, the integral cannot be computed with sufficient accuracy. Using MATLAB's "erf" function (after a coordinate transformation) or a symbolic computation might help.