как переписать код в Matlabe 7

[aniskin]

Ребята у меня проблема с использованием функции fmincon. Раньше в матлабе 6.0 я использовала функцию constr, чтобы минимизировать функционал, зависящий от x0 при ограничениях в виде равенства и неравенства, которые тоже зависят от x0. Ограничения у меня имеют интегральный вид. Подскажите, пожалуйста, могу прислать код.

[photon]

А код сильно большой? Может Вы его тут выставите, используя тег code. Кстати, в версии 6.5 присутствуют обе функции, но в описании constr() сказано, что она замещается функцией fmincon() и в последующих версиях уйдет, поэтому предпочтительнее использовать fmincon(). Поскольку у меня стоит 6.5, то я смогу попробовать поэкспериментировать с обеими функциями - давайте, может что и получится.

[aniskin]

пожалуйста посмотрите, буду очень признательна, если Вы мне поможите.
%****main**********
clear all;
close all;

global n g m lb1 be cc c2 vmax r0 ep dg ZE ZK consta constb constf0 chetchik
pi2=pi*2;
ka=1.4; % dlya vozduha
hg2=(ka+1)/(ka-1);
n=601;
dg=pi2/(n-1);
g=0:dg:pi2;
ep=2;
%beta=10;
be=0.1; %beta*pi/180;
vmax=1;
vinf=1.0;
% lb=0.6
% mb=sqrt(2/(ka+1))*lb/(sqrt(1-(lb^2/hg2)))
mb=0.5;
r0=0.0;
i=sqrt(-1);
pip=pi/2;
ZE=exp(i*g);
ZK=r0*exp(i*(pi+2*be));
%Kochin-Loicans turbulent PS
constf0=-2.00;
consta=1.17;
constb=4.75;

if mb>0;
c2=0.296
lb=sqrt((ka+1)/2)*mb/sqrt(1+(ka-1)*mb*mb/2) %%???
lb1=2*lb/(1+sqrt(1+4*c2*lb*lb));
else
lb=1.0
lb1=1.0
c2=0.0
end
cc=c2*lb1^2/(1+c2*lb1^2);
options(1)=0;
options(2)=1e-4;
options(3)=1e-4;
%options(14)=1e+10;

m=16;
for k=1:m;
x0(k*2-1)=0.0;
x0(k*2)=0.0;
end

x1=fmincon('fun',x0,options);
%************

function [f,f1]=fun(x0)

global n g m lb1 be cc c2 vmax r0 ep dg ZE ZK consta constb constf0 chetchik %constg0 constb %

i=sqrt(-1);
pip=pi/2;

%generate HI,p,q

HI=log(lb1)+((1-ep)-4*i*exp(i*be)*sin(be)*cc-ZK)./ZE-log(1-ZK./ZE);

for k=1:m;
r=(k+1)^(0);
HI=HI+r*(x0(2*k-1)+i*x0(2*k))./(ZE.^(k+1));
end

p=real(HI);
q=imag(HI);

% calculate
ds=exp(-p).*(2*sin((g)/2)).^(ep-1);
ds=ds-c2*exp(p).*(2*sin(g/2)).^(3-ep).*(2*cos(g/2-be)).^2;

dx=ds.*cos(-q+pip+g-(g-pi)*(ep-1)/2);
dy=ds.*sin(-q+pip+g-(g-pi)*(ep-1)/2);

x(1)=0.0;
y(1)=0.0;
s(1)=0.0;

for k=2:length(g),
s(k)=s(k-1)+(ds(k)+ds(k-1)).*(dg/2);
x(k)=x(k-1)+(dx(k)+dx(k-1)).*(dg/2);
y(k)=y(k-1)+(dy(k)+dy(k-1)).*(dg/2);
end;
v=exp(p).*(2*cos(g/2-be)).*(2*sin(g/2)).^(2-ep);
v=v./(1-c2*v.*v);
vm=max(abs(v));

z=x+i*y;
per=max(s);
%**********************
[cmin,imin]=min(abs(v(3:n-1)));
ista=imin;
vv=abs(v);
ssta=s(ista);

form(ista)=0;
for k=ista-1:-1:1, % на верхней стороне
form(k)=form(k+1)+(vv(k)^(constb-1)+vv(k+1)^(constb-1))*abs((s(k)-s(k+1)))/2;
end;

form1(ista)=0;
for k=ista+1:n, % на нижней поверхности
form1(k)=form1(k-1)+abs((s(k)-s(k-1)))*(vv(k)^(constb-1)+vv(k-1)^(constb-1))/2;
end;
% derivativ on upper side ***************
dv(ista)=(v(ista-1)-v(ista))/abs((s(ista-1)-s(ista)));
for k=ista-1:-1:2;
dv(k)=(v(k-1)-v(k+1))/abs(s(k-1)-s(k+1));
end
dv(1)=(v(1)-v(2))/abs(s(1)-s(2));

% derivativ on bottom side*****************
dv1(ista)=(v(ista+1)-v(ista))/abs(s(ista+1)-s(ista));
for k=ista+1:n-1;
dv1(k)=(v(k+1)-v(k-1))/abs(s(k+1)-s(k-1));
end
dv1(n)=(v(n)-v(n-1))/abs(s(n)-s(n-1));

% %****************************************

% form-param. on uper side

form(ista)=consta/constb;
for k=ista-1:-1:1
form(k)=consta*dv(k)*form(k)/(vv(k)^constb);
end
% %****************************************

% form-parametr on bottom side

form1(ista)=consta/constb;
for k=ista+1:n,
form1(k)=consta*dv1(k)*abs(form1(k))/(vv(k)^constb);
end
fp=form1;

%**********************************************

%**********************


f=max(s);
f1(1)=abs(z(n));
f1(2)=abs(vm-vmax);
f1(3)=constf0-min(form(ista:-1:2));
f1(4)=constf0-min(fp)

[aniskin]

.......... ПОГОДИ, КОД НЕМНОГО НЕ ВЕРНЫЙ, НЕ ТО ВЫСТАВИЛА, чуть позже вышлю, спасибо

[photon]

А если Вас устраивала функция constr(), почему Вы не можете просто перетащить ее из старой версии - она не built-in, так что это легко должно быть.

[photon]

А еще в help-e к Matlab 6.5 нашел образец перехода от функции constr() к fmincon(). Приведу этот раздел help-a полностью

Цитата:

Example of Converting from constr to fmincon
Old Call to constr

OPTIONS = foptions;

OPTIONS(13) = 2;             % Two equality constraints
OPTIONS(1) = 1;
OPTIONS(9) = 1;
A1 = [ 1 4 -3]; b1 = 2;
A2 = [ 2 5 0]; b2 = 9;
x0 = [1; .5; .8];
LB = []; UB = [];
[X,OPTIONS,LAMBDA,HESS] = ...
   constr('myfuncon',x0,OPTIONS,LB,UB,'mygradcon',A1,b1,A2,b2);

% myfuncon.m
[F, C] = myfuncon(x,A1,b1,A2,b2)
F = x(1) + 0.0009*x(2)^3 + sin(x(3));

C(1,1) = A1*x-b;             % Equality linear constraint
C(2,1) = 3*x(1)^2-1;         % Equality nonlinear constraint

C(3,1) = A2*x-b2;            % Inequality linear constraint
C(4,1) = 7*sin(x(2))-1;      % Inequality nonlinear constraint

% mygradcon.m
[G, DC] = mygradcon(x,alpha)
G = [1;                       % Gradient of the objective
     3*0.0009*x(2)^2;
     cos(x(3))];

DC(:,1) = A1';                % Gradient of the constraints
DC(:,2) = [6*x(1); 0; 0];
DC(:,3) = A2';
DC(:,4) = [0; 7*cos(x(2)); 0];

Цитата:

New Call to fmincon

OPTIONS = optimset(...
   'Display', 'iter', ...
   'GradCheck', 'on', ... % Check gradients.
   'GradObj', 'on', ...   % Gradient of objective is provided.
   'GradConstr', 'on');   % Gradient of constraints is provided.

A1 = [ 1 4 -3]; b1 = 2;   % Linear equalities
A2 = [ 2 5 0]; b2 = 9;    % Linear inequalities

x0 = [1; .5; .8];
LB = []; UB = [];
[X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = ...
   fmincon(@myfun,x0,A2,b2,A1,b1,LB,UB,@mycon,OPTIONS);

% myfun.m
function [F,G] = myfun(x)
F = x(1) + 0.0009*x(2)^3 + sin(x(3));
G = [1; 
     3*0.0009*x(2)^2;
     cos(x(3))];

% mycon.m
function [C,Ceq,DC,DCeq]= mycon(x)
Ceq(1,1) = 3*x(1)^2-1;         % Equality nonlinear constraint
C(1,1) = 7*sin(x(2))-1;        % Inequality nonlinear constraint
DCeq(:,1) = [6*x(1); 0; 0];    % Gradient of equality 
                               %    nonlinear constraint
DC(:,1) = [0; 7*cos(x(2)); 0]; % Gradient of inequality 
                               %    nonlinear constraint

[aniskin]

спасибо, но это я уже изучила. Оказывается проблема в не написании, а в том, что ограничения сложные слишком... СПАСИБО PHOTON ТЕБЕ ЗА ВНИМАНИЕ.