Optimization- problem based approach- parameter identification using function solve

My aim is obtain value of x such that function expr in the code is minimized. The error I am getting is Objective must be a scalar OptimizationExpression or a struct containing a scalar OptimizationExpression. I have referred this forum and read some relevant posts. However, I am not able to solve the issue and hence this post. If there is a better of doing like not using a for loop to increase efficiency, etc., what I am attempting to do, then please suggest. Here is my code

A = 332.8668;%constant value
B = 128.6184;%constant value
n = 0.4234;%constant value
epsdot = 0.0001;%constant value
s = table2array(ExpDataMatlab);%using import button in GUI I import the data (file attached) and convert it to double.Represents experimental values
sigmaexp = s(:,2);%represents stress observed experimentally
eps = s(:,1);%represents strain
epsdot1 = s(:,3);%represents strain rate
x = optimvar('x',1);
expr = optimexpr;
sigmanum = @(x)(A+B*(eps.^n))*(1+x*log(epsdot1/epsdot));%represents stress calculated numerically where matlab will keep changing x until "expr" is minimized
for i = 1:length(sigmaexp)
    expr = @(x)((1/i)*sum(abs((sigmaexp-sigmanum)/sigmaexp))*100); 
end
JCproblem = optimproblem;
JCproblem.Objective = expr;%this final line is where I get the error. Hence I have given my further code followed by %
%options = optimset('PlotFcns',@optimplotfval);
%init.x = 0.0098 %initial guess of x;
%[sol,fval,exitflag,output] = solve(JCproblem,init,options);

Answer 1

[Walter Roberson]

A = 332.8668;%constant value
B = 128.6184;%constant value
epsdot = 0.0001;%constant value
n = 0.4234;%constant value

All scalars

eps = s(:,1);%represents strain

That's a column vector

epsdot1 = s(:,3);%represents strain rate

Column

sigmanum = @(x)(A+B*(eps.^n))*(1+x*log(epsdot1/epsdot));%represents stress calculated numerically where matlab will keep changing x until "expr" is minimized

eps is column vector, and n is a scalar, so eps.^n is a column vector. B is scalar and scalar times column vector is column vector. A is scalar and scalar plus column vector is column vector. So (A+B*(eps.^n)) is column vector

epsdot1 is a column vector, epsdot is a scalar, so epsdo1/epsdot is a column vector. log of that is a column vector. I think x is intended to be a scalar and scalar times column vector is column vector.

We now have two column vectors * each other. That is a dimension mismatch. P*Q requires that size(P,2) == size(Q,1) but the second dimension of a column vector is 1 and the first dimension of a column vector is not (usually) 1, so you have a problem. You need to transpose one of them if you want (n x 1) * (n x 1)' -> n x n, or you want (n x 1)' * (n x 1) -> 1 x 1. Or you need to switch to the .* operator, which would give you a column vector.

for i = 1:length(sigmaexp)
    expr = @(x)((1/i)*sum(abs((sigmaexp-sigmanum)/sigmaexp))*100); 

sigmanum is a function handle, and you are subtracting the function handle from the column vector. You cannot use a function handle in arithmetic. The only things you can do with a function handle are store them, pass them, or invoke them -- which requires using () with any parameters inside the ()

If you were to use sigmanum(x) then as explored above, you seem to be expecting a vector output from that function. Column vector sigmaexp minus column vector sigmanum is potentially acceptable. But then you have the / operator with a column vector on the right hand side. That is a legal operation, which would give an n x n matrix where n = size(sigmaexp,1) . abs() of that would be n x n, sum() of that would be 1 x n. So expr would be a function handle that returned a 1 x n vector.

Each iteration of for i you overwrite expr with the exact same content. You might as well only do one iteration.

However... your optimization function is required to return a scalar, but this would return a row vector.

[Venkatesh Madhukar Deshpande]

Thank you very much for the quick response @Walter Roberson

The '.' is very important which I did not understand previously. Also I do not need to loop i from 1 to length(sigmexp) since I need sum of all the errors and minimize it later. Here is the updated code

A = 332.8668;
B = 128.6184;
n = 0.4234;
C = 0.0098;
epsdot = 0.0001;
s = table2array(ExpDataMatlab);
sigmaexp = s(:,2);
eps = s(:,1);
epsdot1 = s(:,3);
x = optimvar('x',1,'LowerBound',0.0061,'UpperBound',0.0134);
initialpt.x = 0.0098;
expr = optimexpr;
sigmanum = @(x)(A+B*(eps.^n)).*(1+x*log(epsdot1/epsdot));%first correction made here in the code
i = length(sigmaexp); %2nd change made here
expr = (1/i)*sum((sigmaexp-sigmanum(x))./sigmaexp)*100;%3rd correction made here
JCproblem = optimproblem;
JCproblem.Objective = expr;
%options = optimset('PlotFcns',@optimplotfval);
options = optimoptions(JCproblem,'Display','iter','Algorithm','interior-point');
[sol,fval,exitflag,output] = solve(JCproblem,initialpt,'Options',options);

I have just one question regarding solve. If you want me to put it as separate question, I will do so. I would like to use fmincon to solve my problem i.e. minimize the JCproblem. However, I cannot directly use fmincon instead of solve since it gives me an error saying that JCproblem should be a function (and here it is an optimization expression). Can you suggest how can I invoke fmincon using solve or how can I directly use fmincon function?

[Walter Roberson]

A = 332.8668;
B = 128.6184;
n = 0.4234;
C = 0.0098;
epsdot = 0.0001;
s = table2array(ExpDataMatlab);
sigmaexp = s(:,2);
eps = s(:,1);
epsdot1 = s(:,3);
x = optimvar('x',1,'LowerBound',0.0061,'UpperBound',0.0134);
initialpt.x = 0.0098;
expr = optimexpr;
sigmanum = @(x)(A+B*(eps.^n)).*(1+x*log(epsdot1/epsdot));%first correction made here in the code
i = length(sigmaexp); %2nd change made here
expr = (1/i)*sum((sigmaexp-sigmanum(x))./sigmaexp)*100;%3rd correction made here
JCproblem = optimproblem;
JCproblem.Objective = expr;
%options = optimset('PlotFcns',@optimplotfval);

%options = optimoptions(JCproblem,'Display','iter','Algorithm','interior-point');
options = optimoptions('fmincon','Display','iter','Algorithm','interior-point');


[sol,fval,exitflag,output] = solve(JCproblem, initialpt, 'Options', options, 'solver', 'fmincon');

[Walter Roberson]

Note that the problem you set here can be converted entirely into a linear expression, so there is little benefit to forcing it to use fmincon. linprog can calculate the result quite directly.

[Venkatesh Madhukar Deshpande]

Yes, you are right. It is a linear expression. Actually I want to check if I use any other solver, will it give me the same result. Using linprog the problem does get solved.

It is getting presolved actually. The linear expression expr = -732.432*x + 4.6756. When I use x = 0.0134 (upper bound value), expr = -5.138 and at this point, it says solution is obtained as minimum value is obtained. However, I need to make expr = +732.432*x+4.6756 so that I can minimize this value. Is there anyway I can use absolute value of expr as it puts x values one by one through many iterations and minimize it . In my original code, I used

expr = (1/i)*sum(abs((sigmaexp-sigmanum(x))./sigmaexp))*100;

However, it gave me an error as abs cannot be applied to a function

[Venkatesh Madhukar Deshpande]

I think I have found a workaround to my issue of finding absolute value. Instead of making it positive by using abs, I can square the value and minimize it. expr represents error between experimental and iteratively obtained values and the minimization of the squared expression will also lead to the same x value. So the new expression is

expr = ((1/i)*sum(((sigmaexp-sigmanum(x))./sigmaexp)))^2*(100);
JCproblem = optimproblem;
JCproblem.Objective = expr;
options = optimoptions(JCproblem,'Display','iter');
[sol,fval,exitflag,output] = solve(JCproblem,initialpt,'Options',options);

The expression is a quadratic expression in x and is solved using lsqlin. Thank you for the help @Walter Roberson

[Venkatesh Madhukar Deshpande]

I also changed the subject line of my question so that in future it might be helpful for someone who wants to identify parameter using solve function

[Walter Roberson]

You can use norm() instead of abs()

ExpDataMatlab = table(randi(19,10,1),randi([1 5],10,1),randi([1 9],10,1))

ExpDataMatlab = 10×3 table
    Var1    Var2    Var3
    ____    ____    ____

     12      1       9  
      3      3       7  
     15      3       8  
     19      1       1  
     18      2       8  
      5      5       1  
     11      2       2  
     17      1       5  
     18      2       2  
      6      2       8  

A = 332.8668;
B = 128.6184;
n = 0.4234;
C = 0.0098;
epsdot = 0.0001;
s = table2array(ExpDataMatlab);
sigmaexp = s(:,2);
eps = s(:,1);
epsdot1 = s(:,3);
x = optimvar('x',1,'LowerBound',0.0061,'UpperBound',0.0134);
initialpt.x = 0.0098;
expr = optimexpr;
sigmanum = @(x)(A+B*(eps.^n)).*(1+x*log(epsdot1/epsdot));%first correction made here in the code
i = length(sigmaexp); %2nd change made here

which norm(sigmaexp-sigmanum(x))

/MATLAB/toolbox/optim/problemdef/+optim/+problemdef/@OptimizationExpression/norm.p  % optim.problemdef.OptimizationExpression method

expr = (1/i)*sum(norm(sigmaexp-sigmanum(x))./sigmaexp)*100;%3rd correction made here
JCproblem = optimproblem;
JCproblem.Objective = expr;
%options = optimset('PlotFcns',@optimplotfval);

%options = optimoptions(JCproblem,'Display','iter','Algorithm','interior-point');
options = optimoptions('fmincon','Display','iter','Algorithm','interior-point');


[sol,fval,exitflag,output] = solve(JCproblem, initialpt, 'Options', options, 'solver', 'fmincon');

Solving problem using fmincon.
                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       1    1.426490e+05    0.000e+00    7.017e+05
    1       2    1.376205e+05    0.000e+00    3.711e+05    3.681e-03
    2       3    1.375953e+05    0.000e+00    4.259e-01    1.841e-05

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.

[Venkatesh Madhukar Deshpande]

When I am running code with this expression,

expr = (1/i)*sum(norm(sigmaexp-sigmanum(x))./sigmaexp)*100

it gives me error using norm, first argument must be single or double. I am not able to troubleshoot. sigmaexp is in double form.

[Walter Roberson]

Above that line insert

which norm(sigmaexp-sigmanum(x))

and see what you get out. I get

/MATLAB/toolbox/optim/problemdef/+optim/+problemdef/@OptimizationExpression/norm.p  % optim.problemdef.OptimizationExpression method

at least in R2021a.

You can potentially code

norm(sigmaexp-sigmanum(x))

as

sqrt((sigmaexp-sigmanum(x)).*conj(sigmaexp-sigmanum(x)))

[Venkatesh Madhukar Deshpande]

When I insert

which norm(sigmaexp-sigmanum(x))

Then i get

built-in (C:\Program Files\MATLAB\R2020b\toolbox\matlab\matfun\norm)

I am using Matlab R2020b. With this

sqrt((sigmaexp-sigmanum(x)).*conj(sigmaexp-sigmanum(x)))

it solves the problem. Thank you.