function with myblackbox using fminunc

Hello guys! I got a Function F(y(x)) = sum (( yref-y(x))^2) and x(1) = q and x(2)=r and x=[q;r] and yref=0. I wanted to code this function to be used in a multi-objective optimization etc.

my initial idea is:

function F = myblackbox(x)
q = x(1)
r = x(2);
yref = 0;
y = solvemyoptimizationproblem(q,r);
F = somefunctionofy(y);

but i don't know how to use fminunc here to do a blackbox optimization and how to replace those things to have F(y(x)) correctly.

Answer 1

[Matt J]

lsqnonlin would be better suited to this,

x0=[q_guess,r_guess];

x=lsqnonlin( @(x) yfunction(x(1),x(2))-yref, x0);

[Ali Esmaeilpour]

would you please write complete function or code? because when I want to run your answer it says undefined 'yofx'

[Matt J]

The code is as complete as your description allows. In your problem statement

F(y(x)) = sum (( yref-y(x))^2)

You haven't told us what y(x) looks like.

[Ali Esmaeilpour]

y is closed-loop response and x=[q;r] and q and r are weighted matrices.

[Ali Esmaeilpour]

my closed-loop response contains a 1 by 100 matrix of numerical values and x is a decision variable here and x=[q;r] and I'm going to tune those q and r

[Matt J]

Then yfunction, in my post, should be the function that takes x as input and returns the corresponding 1x100 response.

[Ali Esmaeilpour]

how can I find a function for my closed-loop response y? I got closed-loop response:

y = [-2,-2.06844254709663,-2.10111830954110,-2.10207250524360,-2.07526182090699,-2.02451336009587,-1.95348936029213,-1.86566104239566,-1.76428609501489,-1.65239283602609,-1.53276900719044,-1.40795510853059,-1.28024204134041,-1.15167219215486,-1.02404365938112,-0.898917175767882,-0.777625266347119,-0.661283291728917,-0.550801952158784,-0.446900918848503,-0.350123272116316,-0.260850455112571,-0.179317482440462,-0.105628171524862,-0.0397701935839232,0.0183702310931665,0.0689941424455793,0.112376278061898,0.148852238329798,0.178806457148514,0.202661035909770,0.220865478976649,0.233887353198917,0.242203878667164,0.246294444722945,0.246634033866838,0.243687526483402,0.237904851041483,0.229716937775089,0.219532428513973,0.207735091324267,0.194681885774748,0.180701622852288,0.166094162756435,0.151130093813326,0.136050836531366,0.121069118234489,0.106369765666474,0.0921107653476128,0.0784245442470077,0.0654194263502996,0.0531812239571451,0.0417749259245079,0.0312464485645638,0.0216244182313353,0.0129219582599530,0.00513845621642096,-0.00173870934149472,-0.00773250023679152,-0.0128746228204809,-0.0172040168372426,-0.0207654382076332,-0.0236081426698606,-0.0257846749601092,-0.0273497662907451,-0.0283593408257317,-0.0288696309906560,-0.0289363993725932,-0.0286142641583750,-0.0279561239979665,-0.0270126774854485,-0.0258320316788283,-0.0244593937199454,-0.0229368391498630,-0.0213031503020680,-0.0195937183993203,-0.0178405024207615,-0.0160720382595013,-0.0143134917503068,-0.0125867493914005,-0.0109105408257964,-0.00930058753949017,-0.00776977247844266,-0.00632832582476125,-0.00498402235945812,-0.00374238646270093,-0.00260690103919371,-0.00157921715665090,-0.000659361544098163,0.000154060515616681,0.000863668978079722,0.00147312127494728,0.00198693439542896,0.00241031792821887,0.00274901888833613,0.00300917888601873,0.00319720397995884,0.00331964734521793,0.00338310467267475,0.00339412212095973];

[Ali Esmaeilpour]

I give you the code that I got this close-loop response from and mux1 adn mux2 are closed-loop responses I mean y=mux1 in the following program:

clc;
clear;
close all;
%% Time structure 
t0 = 0;         
tf = 10;        
dt = 0.1;      
time = t0:dt:tf-dt; 
%% System structure  ===> x(k+1) = A*x(k) + B*u(k) + G*w(k)
A = [1.02 -0.1;0.1 0.98];
B = [0.5 0;0.05 0.5];
G = [0.3 0;0 0.3];
[nx,nu] = size(B);
%% MPC parameters
Q = eye(nx);
R = eye(nu)*50;
Qf = eye(nu)*50;
N = 26;
%% Problem parameters
M = blkdiag(Q,R);
w1 = randn(1,100); 
w1 = w1-mean(w1);
muw1 = mean(w1);
sigmaw1 = var(w1);
w2 = randn(1,100); 
w2 = w2-mean(w2);
muw2 = mean(w2);
sigmaw2 = var(w2);
w = [w1;w2];
sigmaw = [sigmaw1 0 ;0 sigmaw2];
x = randn(2,1);
mux{1} = [-2;2];
sigmax{1} = [1,0;0,1];
h1x = [-1/sqrt(5);-2/sqrt(5)];
h2u = [-0.4/sqrt(1.16);1/sqrt(1.16)];
epsilon = 0.5;
g1 = 3;
g2 = 0.2;
c = 0;
%% Main loop
for t = t0:dt:tf-dt
c = c+1;  
sigmax = sdpvar(repmat(nx,1,N+1),repmat(nu,1,N+1));
p = sdpvar(repmat(4,1,N-1),repmat(4,1,N-1));
pN = sdpvar(2,2);
F = sdpvar(repmat(2,1,N+1),repmat(2,1,N+1));
K = sdpvar(repmat(2,1,N+1),repmat(2,1,N+1));
muu = sdpvar(repmat(nu,1,N),repmat(1,1,N));
sigmau = sdpvar(repmat(2,1,N),repmat(2,1,N));
constraint = [];
objective = 0;
for k = 1:N-1   
    mux{k+1} = A*mux{k} + B*muu{k};        
    objective = objective + 1/2*(trace(M*p{k}));         
    constraint2 = [sigmax{k+1} (A*sigmax{k}+B*F{k}) G*sigmaw;(A*sigmax{k}+B*F{k})' sigmax{k} zeros(2);(G*sigmaw)' zeros(2) sigmaw]>=0;
    constraint3 = [sigmau{k} (F{k});(F{k})' sigmax{k}]>=0;
    constraint4 = h1x'*mux{k}<=(1-(0.5*epsilon))*g1-((0.95/2*epsilon*g1*(1-0.95))*h1x'*sigmax{k}*h1x);
    constraint5 = h2u'*muu{k}<=(1-(0.5*epsilon))*g2-((0.95/2*epsilon*g2*(1-0.95))*h2u'*sigmau{k}*h2u);
    constraint6 = [p{k} [sigmax{k};F{k}] [mux{k};muu{k}];[sigmax{k};F{k}]' sigmax{k} zeros(2,1);[mux{k};muu{k}]' zeros(1,2) ones(1)]>=0;
    constraint7 = [F{k} K{k};sigmax{k} eye(2)];
    constraint = [constraint,constraint2,constraint3,constraint4,constraint5,constraint6,constraint7];          
end        
objective = objective + (1/2*trace((trace(Q*pN))));
constraint8 = h1x'*mux{N}<=(1-(0.5*epsilon))*g1-((0.95/2*epsilon*g1*(1-0.95))*h1x'*sigmax{N}*h1x);
constraint9 = [pN-sigmax{N} mux{N};mux{N}' 1];
constraint = [constraint,constraint8,constraint9];
option = sdpsettings('solver','sedumi');
a = optimize (constraint,objective,option);   
u = value(muu{1});
U1(c) = u(1,1);
U2(c) = u(2,1);
U3 = [U1;U2];   
mux3 = value(mux{1});
mux1(c) = mux3(1,1);
mux2(c) = mux3(2,1);
F = value(F{1});
X = value(sigmax{1});
f1{c} = F;
X1{c} = X;
K = F*inv(X);
% K=value(K{1})
K1{c} = K;
K2 = value(K1{c});
x = value(x);
u1 = u+K2*(x-mux3);
u2{c} = u1;
x = A*x+B*u1+G*w(:,c);
mux{1} = value(mux{2});
sigmax{1} = value(sigmax{2});    
ud = value(u2{c});
ud1(c) = ud(1,1);
ud2(c) = ud(2,1);
x1(c) = x(1,1);
x2(c) = x(2,1);       
end
%% Plot results
figure(1)
subplot(4,1,1)
plot(time,ud1,'LineWidth',2)
xlabel('time')
ylabel('ud1')
title('Input signal')
grid on
subplot(4,1,2)
plot(time,ud2,'LineWidth',2)
xlabel('time')
ylabel('ud2')
grid on
subplot(4,1,3)
plot(time,x1,'LineWidth',2)
xlabel('time')
ylabel('x1')
title('State 1')
grid on
subplot(4,1,4)
plot(time,x2,'LineWidth',2)
xlabel('time')
ylabel('x2')
title('State 2')
grid on
figure(2)
subplot(4,1,1)
plot(time,x1,'r',time,mux1,'g','LineWidth',2)
xlabel('time')
legend({'x1','Mean of x1'},'Location','northeast')
grid on
subplot(4,1,2)
plot(time,x2,'r',time,mux2,'g','LineWidth',2)
xlabel('time')
legend({'x2','Mean of x2'},'Location','northeast')
grid on
subplot(4,1,3)
plot(mux1,mux2,'LineWidth',2)
title('State space phase plane')
xlabel('mux1')
ylabel('mux2')
grid on
subplot(4,1,4)
plot(U1,U2,'LineWidth',2)
title('Control space phase plane')
xlabel('U1')
ylabel('U2')
grid on
figure
subplot(4,1,1)
plot(time,mux1,'b','LineWidth',2)
title('Mean of first state')
xlabel('Time')
ylabel('mux1')
grid on
subplot(4,1,2)
plot(time,mux2,'b','LineWidth',2)
title('Mean of second state')
xlabel('Time')
ylabel('mux2')
grid on
subplot(4,1,3)
plot(time,U1,'r','LineWidth',2)
title('First input')
xlabel('Time')
ylabel('U1')
grid on
subplot(4,1,4)
plot(time,U2,'r','LineWidth',2)
title('Second input')
xlabel('Time')
ylabel('U2')
grid on

[Ali Esmaeilpour]

for this code you should have YALMIP in your matlab path and also sedumi as solver installed

[Matt J]

how can I find a function for my closed-loop response y?

To use lsqnonlin or any other Matlab optimizer, all you need is code that generates y for a given choice of [x(1), x(2)], which you seem to have already. For lsqnonlin, you also need this function to be differentiable, and it is not clear from the complexity of your code if that is the case.

SInce we are not sure it is differentiable, then you could use fminsearch instead, which doesn't use derivatives, and works well for a small number of unknowns,

x=fminsearch(@(x) norm(yfunction(x)-yref), [q_guess, r_guess])

[Ali Esmaeilpour]

so I put that fminsearch at the end of my main code?

[Matt J]

fminsearch will search for the optimal x. You put it wherever you need the optimization to occur.